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1  Objectives 

The  general  objectives  of  this  effort,  originally  envisioned  as  a  three-year  project,  are  the 
following; 

1.  To  advance  the  state-of-the-art  of  direct  numerical  simulation  (DNS)  toward  the  ulti¬ 
mate  goal  of  configuration  DNS. 

2.  To  study  the  physics  of  transition  in  hypersonic  wall-bounded  flows.  In  particular,  to 
identify  likely  instability  modes  and  paths  to  transition  and  to  assess  the  effects  on 
transition  of  adverse  pressure  gradient,  wall  cooling,  and  body  asymmetry. 

3.  To  develop  a  transitional-flow  data  base  for  hypersonic  flows  that  can  be  used  1) 
to  extract  transition  physics  through  postprocessing  and  computational  flow  imaging 
(CFI)  and  2)  to  develop  and  refine  transition  models  (not  part  of  this  effort). 

4.  To  begin  development  of  a  configuration  large-eddy  simulation  (LES)  capability  for 
high-speed  transitional  flows  with  benchmark  DNS  calculations  against  which  future 
LES  data  can  be  compared.  (LES  work  per  se  not  currently  part  of  this  effort.) 

It  was  proposed  that  these  general  objectives  be  met  by  considering  two  specific  test 
problems. 

1.  In  the  first  year,  the  laminar-turbulent  transition  of  the  boundary  layer  along  an  ax- 
isymmetric  flared-cone  in  Mach  6  flow  was  considered.  The  intent  was  to  study  the 
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effects  on  transition  of  concave  surface  curvature  and  adverse  streamwise  pressure  gra¬ 
dient. 

2.  The  specific  goal  for  the  two  follow-on  years  was  to  simulate  a  transitional  hypersonic 
boundary-layer  flow  along  on  a  cone  with  an  elliptical  cross-section  at  zero  angle  of 
attack.  Such  a  configuration  is  of  particular  interest  because  it  represents  the  proto¬ 
type  forebody  for  future  hypersonic  flight  vehicles  (e.g.,  the  National  Aerospace  Plane 
(NASP)). 


2  Status  of  Project 

The  specific  first-year  objective  was  met  and  significant  progress  was  made  toward  the  specific 
and  general  objectives  intended  for  the  second  and  third  years  of  the  project. 

Specifically,  progress  to  date  has  proceeded  on  three  fronts,  each  of  which  expands  ex¬ 
isting  DNS  capabilities  in  the  direction  of  “configuration  DNS”  stated  in  general  objective 
1  above.  Previous  experimental  and  numerical  work  on  the  elliptic-cone  problem  by  others 
has  shown  that  crossflow  instabilities,  rather  than  first-  or  second-mode  instabilities,  are 
likely  to  dominate  the  transition  process  when  the  eccentricity  of  the  cross-sectional  ellipse 
is  moderate  to  large.  Thus,  the  problem  is  particularly  difficult  for  two  reasons:  1)  existing 
DNS  capabilities  must  be  expanded  in  terms  of  configuration  (complex  geometry),  and  2) 
little  or  no  experience  exists  in  the  simulation  of  crossflow  instability  for  high-speed  flows. 
Because  of  the  difficulty  of  the  problem,  it  was  deemed  best  to  proceed  incrementally,  from 
the  current  state  of  the  art,  along  three  fronts:  1)  simulation  of  transition  on  a  flared,  ax- 
isymmetric  cone  in  hypersonic  flow,  2)  simulation  of  crossflow  instability  on  a  supersonic 
swept  wing,  and  3)  simulation  of  a  transitional  hypersonic  flow  on  a  cone  with  an  elliptical 
cross  section.  Simulation  1)  was  completed  and  is  thoroughly  document  below.  Problem  2) 
was  addressed,  is  largely  completed,  and  is  also  documented  thoroughly.  Problem  3)  remains 
in  the  formative  stages,  but  important  groundwork  was  laid.  The  specific  accomplishments 
for  each  of  these  problems  are  discussed  below. 


3  Accomplishments 

3.1  DNS  of  Instability  and  Transition  on  a  Flared  Cone 

INTRODUCTION:  The  computation  documented  in  the  recent  papers  of  Pruett  et  al. 
[15]  (Attachment  1),  Pruett  and  Chang  [16],  and  Pruett  and  Zang  [17],  represents  the  first  at¬ 
tempt  to  numerically  simulate  the  transition  of  a  high-speed  boundary-layer  flow  for  Reynolds 
numbers  and  configurations  of  engineering  interest.  Specifically,  their  work  examined  the 
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laminar  breakdown  of  a  Mach  8  boundary-layer  flow  on  an  axisymmetric  sharp  cone  in  the 
absence  of  streamwise  pressure  gradient.  The  computation  approximately  simulated  the 
sharp-cone  experiment  of  Stetson  et  al.  [21]  and  was  motivated  by  the  hope  that,  when 
experiment  and  computation  are  both  brought  to  bear  on  the  same  physical  problem,  the 
results  will  be  more  enlightening  than  when  either  approach  alone  is  used.  In  the  present 
paper,  we  extend  the  capability  for  “configuration  DNS”  to  a  flared-cone  configuration. 

Recently,  Wilkinson  and  Anders  [24],  Lachowicz  et  al.  [9],  and  Kimmel  [8]  have  each 
conducted  hypersonic  wind-tunnel  experiments  on  axisymmetric  flared-cone  models  to  assess 
the  direct  and  indirect  effects  on  transition  of  streamwise  surface  curvature.  We  consider 
here  only  the  case  of  concave  curvature.  A  concave  flare  induces  several  effects:  an  adverse 
(positive)  streamwise  pressure  gradient,  a  decrease  in  mean  edge  Mach  number,  an  increase 
in  mean  wall  temperature,  a  reversal  of  boundary-layer  growth,  and  Goertler  instability.  The 
conventional  wisdom,  based  largely  on  experience  with  low-speed  flows,  is  that  the  adverse 
streamwise  pressure  gradient  should  promote  more  rapid  transition  than  for  a  configuration 
without  a  flare.  The  results  of  Kimmel  [8]  are  inconclusive,  however.  In  high-speed  boundary 
layers,  the  transition  zone  can  be  quite  long,  due  to  the  generally  stabilizing  influence  of 
increasing  Mach  number.  Whereas  adverse  pressure  gradient  appears  to  accelerate  transition 
onset,  in  some  cases  the  transition  zone  was  observed  to  lengthen  relative  to  the  zero  pressure 
gradient  case.  Apparently,  at  high  speeds,  there  is  interplay  among  subtle  effects,  several  of 
which  were  mentioned  above.  For  example,  possible  interaction/ competition  between  second¬ 
mode  [12]  and  Goertler  instabilities  may  occur,  the  latter  of  which  are  linearly  stable  for  cones 
at  zero  angle  of  attack  without  flare.  Moreover,  variations  in  edge  Mach  number  and  wall 
temperature  may  tend  to  select  or  deselect  certain  modes  of  instability  in  different  regions  of 
the  transition  zone.  A  further  subtlety  concerns  the  level  of  freestream  turbulence  and  the 
issue  of  boundary-layer  receptivity,  which  is  beyond  the  scope  of  the  current  investigation. 
However,  one  should  keep  in  mind  that,  whereas  the  experiment  of  Lachowicz  et  al.  [9]  was 
conducted  in  a  low- turbulence  “quiet”  facility,  the  experiments  of  Kimmel  [8]  were  conducted 
in  a  conventional  facility  with  a  relatively  high  level  of  freestream  turbulence.  Suffice  it  to 
say,  transitional  high-speed  flows  remain  poorly  understood,  and  further  experiment  and 
computation  are  warranted. 

In  the  present  work,  we  adapt  the  approach  and  algorithm  of  Pruett  et  al.  [15]  to  simulate 
the  early  stages  of  laminar-turbulent  transition  on  an  axisymmetric  flared  cone  in  Mach  6 
flow.  Our  motivations  are  to  advance  the  state-of-the-art  in  configuration  DNS  and  to 
study  the  subtle  transition  phenomena  associated  with  surface  curvature  mentioned  above. 
Consequently,  the  effects  of  streamwise  surface  curvature,  mean  streamwise  pressure  gradient, 
and  varying  boundary-layer  edge  conditions  are  incorporated  into  the  governing  equations 
and  boundary  conditions.  A  DNS  test  case  has  been  designed  based  on  the  configuration 
and  flow  conditions  of  the  experimental  investigation  of  Lachowicz  et  al.  [9].  This  choice 
of  test  case  affords  the  additional  opportunity  to  draw  upon  the  theoretical  investigation  by 
Balakumar  and  Malik  [1],  who  conducted  N-factor  studies  for  the  flared-cone  geometry  of 
the  Lachowicz  et  al.  experiment. 
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In  the  next  sections,  we  discuss  the  governing  equations,  present  the  test  case,  summarize 
the  numerical  methodology,  present  results,  and  draw  some  conclusions,  respectively. 

GOVERNING  EQUATIONS:  Let  p,  p,T,  p,  and  k  denote  the  pressure,  density,  temper¬ 
ature,  viscosity,  and  thermal  conductivity  of  a  compressible  fluid,  respectively.  We  assume 
that  the  flow  of  a  compressible  fluid  is  exactly  governed  by  the  compressible  Navier- Stokes 
equations  (CNSE),  which  can  be  found  in  many  references  on  fluid  mechanics.  If  we  restrict 
attention  to  axisymmetric  bodies  and  define  x  as  the  surface  arc  length,  0  as  the  azimuthal 
angle,  and  2;  as  the  coordinate  normal  to  the  wall,  then  in  body-fitted  coordinates,  the  CNSE 
assume  the  form  given  in  Pruett  et  al.  [15],  where  u,  v,  and  w  denote  the  velocity  compo¬ 
nents  in  the  principal  coordinate  directions,  respectively.  We  note  that  the  equations  given 
in  Pruett  et  al.  include  both  transverse  and  streamwise  curvature  terms  and  incorporate  the 
notation  R{x)  and  4>{x)  to  denote  the  body  radius  and  the  angle  of  the  surface  tangent  rela¬ 
tive  to  the  axis,  repectively,  as  functions  of  arc  length.  Four  dimensionless  parameters  arise 
from  nondimensionalization  of  the  governing  equations:  Mach  number  M,  Reynolds  number 
Re,  Prandtl  number  Pr,  and  the  ratio  of  specific  heats,  7.  Throughout  this  work,  the  Pr 
and  7  are  assumed  to  be  constants  with  values  of  0.72  and  1.4,  respectively.  Finally,  the  the 
dependency  of  viscosity  p  on  temperature  is  modeled  by  Sutherland’s  law,  and  k  =  pjPr. 

TEST  CASE:  The  experiment  of  Lachowicz  et  al.  [9],  on  which  the  numerical  computation 
is  based,  examines  Mach  6  flow  on  a  20-inch-long  flared-cone  model  at  zero  angle  of  attack. 
The  model  consists  of  a  5-degree  half-angle  axisymmetric  sharp  cone  blended  to  an  afterbody 
with  a  circular-arc  flare  beginning  10  axial  inches  from  the  apex.  The  radius  of  curvature  of 
the  flare  is  constant  and  equal  to  93.07  inches;  a  curvature  discontinuity  exists  at  the  flare 
interface.  The  nominal  pre-shock  freestream  conditions  are 

Moo  =  6.0 

=  98.78°  R 

(r;  =  810.0°  R)  (1) 

=  11.86  psf 
{p;  =  18, 720  psf) 

Rel  =  2,728,087  ft~^ 

where  Rti  is  the  freestream  unit  Reynolds  number.  Throughout  this  work,  asterisks  denote 
dimensional  quantities,  and  subscripts  00,  t,  and  e  denote  freestream,  total,  and  boundary- 
layer  edge  conditions,  respectively. 

A  precise  thermal  boundary  condition  at  the  model  surface  cannot  be  ascertained  on 
the  basis  of  information  reported  about  the  experiment.  Run  times  were  long,  and  it  is 
believed  that  the  wall  behaved  approximately  adiabatically  under  equilibrium  conditions. 
However,  the  model  was  thin-skinned  and  uninsulated,  and  mounting  hardware  prevented 
the  attainment  of  thermal  equilibrium  near  the  models’  base.  Given  the  uncertainty  as  to 
the  experimental  thermal  condition  and  the  fact  that  we  consider  only  the  early  stages  of 
transition  for  the  DNS,  we  have  adopted  the  boundary  condition  of  Pruett  and  Chang  [16]. 
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The  use  of  this  thermal  condition  is  advantageous  in  that  it  allows  comparison  of  simu¬ 
lated  flared-cone  and  straight-cone  results  for  the  same  boundary  conditions.  The  boundary 
conditions  will  be  specified  precisely  later. 

METHODOLOGY:  Roughly  speaking,  spatial  DNS  is  a  three-step  process.  A  laminar 
and  steady  base  state  is  defined,  whose  stability  is  to  be  investigated.  The  base  state  is 
subject  to  time-dependent  fluctuations  at  or  near  the  computational  inflow  boundary,  and 
the  spatial  evolution  of  these  disturbances  is  calculated  by  solution  of  the  three-dimensional 
and  unsteady  CNSE.  The  temporally  and  spatially  evolving  numerical  solution  is  sampled 
and  stored  periodically  in  time  for  the  subsequent  postprocessing  and  analysis  that  constitute 
the  final  step.  We  address  each  of  these  steps  in  turn. 

To  obtain  the  laminar  base  state  for  their  stability  calculations,  Balakumar  and  Malik 
[1]  computed  the  flow  over  the  flared-cone  configuration  by  both  Euler  and  Navier- Stokes 
techniques.  The  Navier-Stokes  solution  was  sensitive  to  the  choice  of  limiter  (dissipation) 
and  required  very  high  resolution  within  the  wall  layer  to  provide  accurate  first  and  second 
wall-normal  derivatives,  which  were  needed  for  their  linear-stability  analysis.  They  found 
that  the  pressure  distributions  computed  by  the  two  methods  agreed  very  well  except  near 
the  tip  of  the  cone.  Moreover,  N-factors  computed  by  the  two  methods  agreed  to  within  ap¬ 
proximately  ten  percent.  Because  the  difference  in  N-factors  was  small,  in  the  final  analysis, 
for  computational  efficiency,  they  recommended  an  Euler /boundary-layer  approach  over  the 
solution  of  the  full  CNSE. 

For  the  present  work,  we  also  favor  a  boundary-layer  equation  approach.  Specifically, 
we  exploit  the  spectrally  accurate  boundary-layer  code  of  Pruett  and  Streett  [13]  with  the 
modification  for  varying  edge  conditions  proposed  by  Pruett  [14].  The  post-shock  edge 
pressure  distribution  p*{x*)  needed  by  the  boundary-layer  code  is  obtained  from  the  Navier- 
Stokes  solution  of  Balakumar  and  Malik  [1].  In  the  vicinity  of  the  tip  (a  region  excluded  in 
the  DNS  calculation),  the  pressure  distribution  is  modified  by  extrapolating  linearly  from 
downstream  values.  For  self  consistency,  all  other  boundary- layer  edge  conditions  are  inferred 
from  the  pressure,  based  on  the  assumption  of  isentropic  post-shock  flow.  Although  this 
assumption  is  not  strictly  valid  in  the  flared  region,  the  shock  is  highly  oblique  and  weak, 
and  errors  committed  by  assuming  isentropic  post-shock  flow  are  small.  Figure  1  presents 
the  radius  R*  of  the  cone  as  a  function  of  surface  arc  length  x*.  The  edge  Mach  number, 
edge  temperature,  and  edge  pressure  distributions  are  also  provided  as  functions  of  x*.  Note 
that  the  edge  Mach  number  decreases  almost  linearly  from  nearly  5.6  at  the  beginning  of 
the  flare  to  4.9  at  the  end  of  the  model.  In  contrast,  edge  pressure  and  temperature  both 
increase  significantly  and  nearly  linearly  in  the  flare  region.  The  evolution  of  the  laminar 
boundary-layer  displacement  thickness  is  shown  in  Fig.  2.  Note  the  reversal  of  the  normal 
thickening  of  the  boundary  layer  in  the  flare  region.  Typical  temperature  and  streamwise 
velocity  profiles  for  an  adiabatic  wall  are  shown  in  Fig.  3.  The  profiles  correspond  to  axial 
station  1.1725  ft.  in  the  flare  region.  Note  that  the  velocity  distribution  is  nearly  linear 
across  the  boundary  layer,  a  characteristic  of  high-speed  flow,  and  that,  because  of  frictional 
heating,  the  wall  temperature  is  nearly  six  times  that  at  the  boundary  layer  edge. 
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The  boundary- layer  solution  is  interpolated  onto  a  DNS  grid  using  a  spectrally  accurate 
interpolation  code,  which  extrapolates  outside  the  boundary- layer  region  by  solving  the 
continuity  equation  exactly  in  the  asymptotic  limit  as  z  — >■  oo.  The  asymptotic  outer  solution 
is  matched  to  the  inner  solution  at  the  boundary-layer  edge.  Special  attention  is  paid  to 
accurately  computing  the  wall- normal  velocity,  for  which  we  use  the  method  of  Pruett  [14]. 
A  typical  wall-normal  velocity  distribution  versus  2  is  shown  in  Fig.  4.  For  axisymmetric 
cones,  the  wall-normal  velocity  is  negative,  which  implies  that  the  far-field  computational 
boundary  should  be  treated  as  an  inflow  boundary. 

The  special  care  used  to  match  inner  and  outer  solutions  consistently  is  necessary  to 
prevent  discontinuities  in  higher  derivatives  of  the  basic  state,  which  would  contaminate  the 
subsequent  stability  analysis.  As  a  consistency  check,  we  typically  compute  the  residuals 
of  the  time-independent  and  dimensionless  CNSE  by  using  the  interpolated  boundary-layer 
solution  as  an  initial  state.  Residuals  that  are  quite  small,  say  of  order  10“^,  as  shown 
in  Fig.  5,  indicate  that  the  boundary- layer  solution  closely  approximates  the  exact  solution. 
Because  the  residuals  reflect  the  sum  of  approximation  errors,  interpolation  errors,  numerical 
errors,  and  blunders,  small  residual  errors  provide  a  stringent  test  of  self-consistency  among 
the  various  numerical  tools  exploited  by  DNS.  Indeed,  on  several  occasions  we  have  uncovered 
and  corrected  algorithmic  bugs  simply  by  examining  the  steady-state  residuals. 

For  the  DNS  calculation,  we  adapt  the  techniques  of  Pruett  et  al.  [15].  Their  fully  explicit 
collocation  scheme  exploits  the  third-order  low-storage  Runge-Kutta  scheme  of  Williamson 
[25]  for  time  advancement  and  both  spectral  and  central  compact-difference  methods  for 
spatial  discretization.  Specifically,  spectral-collocation  methods  are  used  in  the  azimuthal 
direction,  in  which  the  flow  is  periodic,  whereas  the  high-order  compact- difference  methods 
of  Lele  [10]  are  exploited  in  the  aperiodic  streamwise  and  wall-normal  directions.  Modest 
filtering  is  used  to  correct  the  innate  tendency  of  central-difference  approximations  toward 
odd-even  decoupling  of  the  solution  in  the  far  field  whenever  resolution  is  marginal.  Filtering 
also  alleviates  numerical  reflections  at  the  outflow  boundary.  A  discussion  of  the  necessity 
for  and  implementation  of  filtering  can  be  found  in  Pruett  et  al.  [15].  With  regard  to 
the  spectral-collocation  scheme  used  in  the  azimuthal  direction,  we  implement  dealiasing 
according  to  the  1  /2-rule  to  suppress  high-frequency  oscillations  that  could  lead  to  numerical 
instability. 

Minor  modifications  to  the  scheme  of  Pruett  et  al.  are  necessitated  by  particulars  of  the 
current  fiared-cone  problem.  Because  boundary-layer  edge  conditions  vary  for  the  current 
problem,  flow  quantities  are  nondimensionalized  by  the  pre-shock  freestream  temperature, 
velocity,  and  density,  etc.,  given  or  implied  in  Eq.  1.  Previously,  flow  quantities  were  scaled 
by  appropriate  post-shock  values  at  the  boundary-layer  edge,  which  were  presumed  constant 
along  the  axisymmetric  cone  without  flare.  For  both  the  work  of  Pruett  et  al.  and  the 
current  work,  lengths  are  scaled  by  the  boundary-layer  displacement  thickness  at  the 
computational  inflow  boundary  Previously,  R{x)  was  computed  internally  and  ^(x)  was 
constant;  for  the  present  problem,  both  quantities  are  computed  externally  by  the  boundary- 
layer  code  and  read  as  input  to  the  DNS  algorithm.  As  was  done  in  Pruett  et  al.,  grid 
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spacing  is  uniform  in  the  x  direction,  and,  in  the  wall-normal  direction,  grid  points  are 
tightly  clustered  near  the  wall  and  near  the  critical  layer.  For  the  present  problem,  the 
DNS  computation  includes  only  a  region  of  flow  aft  of  the  flare,  where  the  boundary-layer 
thickness  diminishes  (Refer  to  Fig.  2).  This  is  in  contrast  to  the  growing  boundary  layer 
along  the  cone  without  flare.  In  either  case,  to  account  for  the  variation  in  boundary-layer 
thickness,  we  exploit  a  mapping  of  the  form 

z  =  Tjf{x)  (2) 

However,  for  the  present  problem,  the  mapping  function  is  linear  and  is  given  by 

f(x)  =  1  -  .05— (3) 
^in  ^out 

where  is  the  coordinate  of  the  computational  outflow  boundary. 

The  boundary  conditions  for  the  present  problem  are  taken  without  modification  from 
Pruett  et  al.  [15].  In  brief,  all  flow  variables  are  specified  along  the  inflow  boundary,  a  buffer 
domain  approach  (Streett  and  Macaraeg  [22])  is  exploited  at  the  outflow  boundary,  and 
Thompson  [23]  conditions  are  adapted  along  the  far-field  boundary  ij  =  T/max-  At  the  wall, 
all  velocity  components  vanish;  the  thermal  boundary  condition  may  be  either  adiabatic  or 
isothermal.  As  mentioned  previously,  to  facilitate  comparison  with  the  straight-cone  results 
of  Pruett  and  Chang  [16],  we  adopt  their  hybrid  thermal  boundary  condition  whereby  wall 
temperature  is  held  fixed  at  its  laminar  adiabatic-wall  value.  Density  at  the  wall  is  computed 
directly  from  the  governing  equations. 

PRIMARY  INSTABILITY:  The  hot-wire  data  of  Lachowicz  et  al.  [9]  were  examined 
to  identify  unstable  disturbances.  These  data  show  that  the  most  unstable  frequency  F* 
downstream  of  the  flare  is  in  the  range  218-228  kHz,  in  close  agreement  with  the  230kHz 
predicted  by  Balakumar  and  Malik  [1]  on  the  basis  N-factor  studies.  Unfortunately,  the 
signal-to-noise  levels  of  the  data  of  Lachowicz  et  al.  are  too  low  in  the  region  ahead  of  the 
flare  to  provide  useful  information  there.  Neither  was  information  available  in  regard  to 
the  obliqueness  angle  of  the  primary  instability.  For  this  work,  additional  information  was 
needed  to  identify  a  transition  mechanism.  Consequently,  a  parameter  study  was  conducted 
by  Chang  [6]  on  the  basis  of  nonlinear  parabolized  stability  equation  (PSE)  methodology  to 
identify  a  likely  transition  scenario.  Chang  considered  primary  disturbances  with  temporal 
circular  frequencies  u*  =  2wF*  close  to  those  observed  or  predicted  in  the  investigations 
of  [9]  or  [1].  The  study  confirmed,  as  anticipated,  that  Goertler  modes  of  instability  are 
linearly  unstable  for  a  broad  range  of  azimuthal  wavenumbers  /5*.  In  contrast,  Goertler 
modes  (which  have  zero  temporal  frequency)  are  linearly  stable  for  the  cone  without  flare 
and  gain  energy  only  through  nonlinear  interactions.  F^or  F*  —  230  kHz,  a  fundamental 
azimuthal  wavenumber  given  by  n  =  /3*R*  =  15  was  shown  to  lead  to  highly  unstable 
disturbances.  (Note  that  because  of  periodicity  in  the  azimuthal  direction,  n  must  be  an 
integer.)  A  transition  scenario  was  constructed  based  on  a  disturbance  consisting  initially 
of  three  fundamental  components  whose  nonlinear  interactions  lead  to  a  transitional  state. 
Subsequently,  we  refer  to  the  fundamental  components  and  their  harmonics  by  integer  pairs 
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where  the  integers  I  and  m  identify  the  harmonic  with  respect  to  the  fundamental 
values  of  temporal  frequency  u;  and  azimuthal  wavenumber  /3,  respectively. 

The  amplitudes  of  selected  harmonics  from  the  PSE  calculation  are  shown  in  Fig.  6.  The 
PSE  computation  spanned  the  region  x*  =  0.58  ft.,  ahead  of  the  flare,  to  x*  =  1.4  ft.  in  the 
flared  zone.  At  the  inflow  boundary,  the  calculation  was  initialized  with  a  disturbance  triad 
comprised  of  Fourier  components  (1,0),  (1,1),  and  (1,-1)  (not  shown).  The  (1,0)  component 
represents  an  axisymmetric  disturbance,  wherccis  the  (1,1)  and  (1,-1)  components  represent 
symmetric  oblique  (helical)  disturbances.  All  other  harmonics  emerged  spontaneously  from 
the  nonlinear  interactions  among  the  fundamental  triad.  Fig.  6  shows  saturation  of  the 
fundamental  components  of  the  disturbance  near  x*  =  1.0  ft.  occuring  simultaneously 
with  significant  distortion  of  the  mean  flow  (0,0)  and  rapid  growth  of  higher  harmonics 
to  significant  amplitudes.  These  conditions  indicate  that  the  flow  is  transitioning  from  a 
laminar  to  a  turbulent  state. 

RESULTS:  The  computational  domain  of  the  DNS  encompassed  a  section  along  the  body 
downstream  of  the  beginning  of  the  flare.  Specifically,  the  computational  domain  spanned 
the  region 


x;„  =  0.9<  X-  <1.2  =  x-„,ft. 

O'jr 

0  <  e  <—  (n  =  15)  (4) 

n 

0  <  ri  <7.5 

The  time-periodic  inflow  condition  at  Xj^^  =  0.9  was  derived  from  the  PSE  results  at  the  cor¬ 
responding  station.  At  this  station,  the  fundamental  (1,0)  component  of  the  disturbance  had 
attained  a  large  (nonlinear)  amplitude,  and  numerous  harmonics  of  the  disturbance  showed 
significant  amplitudes.  All  significant  harmonics  were  incorporated  into  the  DNS  inflow 
condition.  Relative  to  the  wavelength  of  the  fundamental  components  of  the  disturbance, 
the  computational  domain  encompassed  approximately  26  wavelengths  in  streamwise  extent. 
The  resolution  of  the  computational  grid  was  2304  x  18  x  128  in  the  streamwise,  azimuthal, 
and  wall-normal  directions,  respectively.  At  this  streamwise  grid  density,  each  fundamental 
wavelength  was  resolved  by  approximately  90  grid  points,  sufficient,  according  to  Pruett  et 
al.  [15],  to  resolve  at  least  the  first  six  streamwise  harmonics.  The  azimuthal  resolution, 
significantly  less  than  that  of  Pruett  et  al.,  was  sufficient  only  for  the  initial  stages  of  laminar 
breakdown,  and  was  dictated  largely  by  available  computational  resources.  Consequently, 
the  simulation  was  halted  prior  to  end-stage  transition. 

The  DNS  computation  was  performed  on  the  Waterways  Experiment  Station  (WES) 
Cray  C90  at  Vicksburg,  MS.  Approximately  250  hours  of  single  processor  time  and  115 
megawords  of  main  memory  were  required  for  the  computation.  The  code  is  optimized  to 
rim  in  parallel,  and  during  the  production  runs,  we  typically  requested  three  CPUs.  In 
physical  terms,  the  computation  proceeded  for  approximately  30  periods  of  oscillation  at 
the  fundamental  frequency.  The  dominant  instability  waves  were  observed  to  propagate 
downstream  at  approximately  90  percent  of  the  boundary-layer  edge  velocity;  consequently. 
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the  physical  computation  time  was  approximately  that  necessary  for  the  a  fundamental 
instability  wave  to  complete  one  traverse  of  the  computational  domain. 


During  the  computation,  the  numerical  solution  was  periodically  written  to  output  and 
archived  at  a  rate  of  64  samples  per  period  of  the  fundamental  frequency.  These  data  were 
later  postprocessed  to  extract  quantitative  and  qualitative  information  about  the  evolution 
of  the  flow.  Most  of  the  results  to  follow  were  obtained  by  Fourier  time-series  analysis  of  the 
final  period  (29-30)  of  the  computation. 

Figure  7,  obtained  by  Fourier  transforms  in  the  periodic  temporal  and  azimuthal  dimen¬ 
sions,  displays  the  streamwise  evolution  of  the  amplitudes  of  selected  harmonics  of  the  DNS 
calculation.  The  corresponding  PSE  results  are  also  provided  for  comparison.  Several  impor¬ 
tant  observations  can  be  made.  First,  the  DNS  and  PSE  results  agree  well  for  approximately 
the  first  third  of  the  computational  domain  of  the  DNS.  Downstream  of  the  first  third  of  the 
domain,  the  DNS  and  PSE  results  gradually  diverge,  but  the  methods  continue  to  indicate 
agreement  in  the  qualitative  behavior  of  the  selected  harmonics.  Downstream  of  x*  =  1.15 
ft.,  some  vestige  of  the  initially  undisturbed  laminar  state  remains;  here  the  DNS  results 
cannot  be  considered  fully  developed  and  should  be  disregarded. 


Figure  8  displays  the  streamwise  evolution  of  various  mean  quantities  computed  by  av¬ 
eraging  both  temporally  and  spatially  over  the  azimuthal  coordinate.  For  comparison,  the 
corresponding  quantities  are  shown  also  for  the  undisturbed  initial  laminar  state.  Clear  in¬ 
dications  of  transitioning  flow  are  evident:  principally  a  dramatic  drop  in  the  shape  factor 
H  =  below  its  laminar  value.  We  note  that  the  laminar  shape  factor,  which  is  nearly 

constant  for  the  cone  without  flare  (Pruett  and  Zang  [17]),  diminishes  with  x  in  the  flared 
region  for  the  present  configuration  because  of  the  decrease  in  Mach  number.  At  this  stage 
of  the  transition  process,  the  decrease  in  H  is  due  almost  entirely  to  an  increase  in  momen¬ 
tum  thickness  82  rather  than  to  a  decrease  in  displacement  thickness  8\.  Near  the  end  of 
the  computational  domain,  the  skin-friction  coefficient  C/  has  increased  approximately  25 
percent  over  its  laminar  value,  where 


^wall 


^^wall 

_du, 

'‘a?'-" 


(5) 


However,  at  the  same  station,  the  coefficient  Ch  of  thermal  stress  at  the  wall  has  more  than 
doubled.  Taken  together,  the  plots  for  Cf  and  Ch  call  attention  to  an  ambiguity  in  defining 
the  location  x*^  of  transition  onset  for  high-speed  flows.  On  the  basis  of  mimimum  wall 
shear,  x*^  =  0.95  ft.  On  the  basis  of  minimum  heat  transfer,  transition  onset  apparently 
occurs  upstream  of  the  computational  inflow  boundary. 


An  overall  impression  is  that,  relative  to  the  straight-cone  case,  transition  on  the  flared 
cone,  as  indicated  by  skin-friction  and  heat  transfer,  begins  earlier  and  proceeds  more  grad¬ 
ually.  The  reason  for  this  qualitative  behavior  may  have  to  do  with  the  linear  instability  of 
Goertler  modes  for  the  flared-cone  problem.  These  modes  manifest  significant  velocities  near 
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the  wall,  and  because  of  their  linear  instability,  they  grow  to  significant  amplitudes  early  in 
the  transition  process.  As  a  result,  mean  skin  friction  is  also  affected  relatively  early.  In  con¬ 
trast,  for  the  straight  cone,  the  disturbances  tend  to  “live”  far  from  the  wall  near  the  critical 
layer.  Disturbance  activity  spreads  toward  the  wall  only  through  strongly  nonlinear  effects, 
which  occur  much  later  in  the  transition  process.  This  observation  of  early  transition  onset 
but  gradual  transition  for  the  cone  with  concave  flare  is  not  inconsistent  with  the  experi¬ 
mental  results  of  Kimmel  [8].  However,  because  different  researchers  use  different  measures 
of  transition  onset,  it  is  difficult  to  draw  definitive  comparisons  between  computation  and 
experiment.  In  our  judgement,  the  transition  community  would  benefit  by  consensus  in  the 
definition  of  transition  onset  for  high-speed  boundary-layer  flows. 

Figure  9  presents  the  streamwise  evolution  of  the  principal  components  of  the  Reynolds 
stress  tensor,  namely  pu'u\  pv'v',  and  pw'w',  where  overlines  denote  temporal  and  azimuthal 
averages  and  primes  denote  fluctations  about  the  mean  states.  Dramatic  differences  relative 
to  the  results  of  Pruett  and  Chang  [16]  and  Pruett  and  Zang  [17]  are  immediately  obvious. 
First,  whereas  for  the  cone  without  flaxe,  the  Reynolds  stresses  are  dominated  by  the  stream- 
wise  component,  along  the  flared  cone  the  streamwise  and  wall-normal  stresses  are  of  nearly 
equal  magnitudes.  Second,  the  dominant  Reynolds  stresses  show  profiles  which  peak  quite 
close  to  the  wall,  in  contrast  to  the  peaks  that  occur  about  a  displacement  thickness  from  the 
wall  in  Pruett  and  Chang  [16].  Similar  trends  are  observed  in  turbulent  Mach  number  Mt 
and  (density-weighted)  turbulent  kinetic  energy  pK  as  shown  in  Fig.  10,  each  of  which  peaJrs 
close  to  the  wall  for  the  flared-cone  case.  (Turbulent  kinetic  energy  is  defined  as  one-half  the 
trace  of  the  Reynolds  stress  tensor.) 

The  differences  noted  above  between  the  flared-cone  and  straight-cone  cases  may  be 
due  to  a  combination  of  effects  manifest  in  other  quantities.  Figure  10  also  displays  the 
streamwise  evolution  of  the  root-mean-square  density  fluctuation.  For  the  flared  cone,  the 
density  fluctuation  is  very  large  at  the  wall,  nearly  30  percent  of  the  reference  density.  For 
the  case  of  the  cone  without  flare,  the  density  fluctuations  are  non-zero  but  relatively  small 
at  the  wall.  Both  configurations  show  peaks  near  the  critical  layer  (about  one  displacement 
thickness  from  the  wall).  We  speculate  that  the  large  density  fluctuation  near  the  wall  is 
directly  related  to  the  significant  streamwise  pressure  gradient  of  the  flared-cone  problem  and 
its  associated  in  wall-temperature  gradient  (Fig.  12).  As  shown  in  Fig.  11,  most  harmonics 
of  the  disturbance  manifest  significant  velocity  fluctuations  in  close  proximity  to  the  wall, 
in  contrast  to  the  usual  character  of  second-mode  disturbances.  The  combination  of  large 
density  fluctuations  and  significant  velocity  fluctuations  near  the  wall  probably  accounts 
for  the  trends  observed  above  for  Reynolds  stress  and  (density- weighted)  turbulent  kinetic 
energy. 

In  addition  to  displaying  the  structure  (amplitude  envelope)  of  the  velocity  and  temper¬ 
ature  components  of  selected  harmonics.  Figures  11  and  13  also  present  a  comparison  of 
DNS  and  PSE  results.  The  comparison  is  made  at  x*  =  1.05  ft.,  the  midpoint  of  the  com¬ 
putational  domain  of  the  DNS.  In  contrast  to  the  previous  close  agreement  between  DNS 
and  PSE  obtained  by  Pruett  et  al.  [15]  and  Pruett  and  Chang  [16],  the  agreement  between 
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methods  can  only  be  considered  good  in  a  qualitative  sense  for  the  flared  cone,  except  for  the 
fundamental  (1,0)  component,  for  which  the  agreement  is  relatively  good.  The  discrepancy 
comes  as  a  surprise,  and,  at  the  present,  we  do  not  know  the  cause. 

Finally,  we  note  that  Lachowicz  et  al.  [9]  observes  “second-mode”  waves  at  the  edge 
of  the  boundary  layer  in  schlieren  images  made  during  his  flared-cone  experiment.  Similar 
structures  have  been  observed  in  hypersonic  boundary-layer  flows  by  several  other  experi¬ 
mentalists  and  have  been  referred  to  previously  as  “rope-like”  waves.  See,  for  example,  the 
excellent  review  by  Smith  [20].  Computational  flow  imaging  (CFI)  of  the  present  numeri¬ 
cal  experiment  also  clearly  reveals  instability  waves  of  rope-like  appearance  (Fig.  14).  To 
crudely  approximate  the  physics  involved  in  schlieren  imagery,  the  instantaneous  wall-normal 
density  gradient  has  been  averaged  over  the  azimuthal  coordinate  and  projected  as  a  two- 
dimensional  image  in  greyscale.  The  numerically  derived  image  is  qualitatively  similar  in 
appearance  to  the  actual  schlieren  image  obtained  by  Lachowicz  et  al.  [9].  In  particular,  the 
wavelength  of  the  instability  is  approximately  two  boundary-layer  displacement  thicknesses, 
in  keeping  with  the  nature  of  second-mode  disturbances. 

CONCLUSIONS:  The  early  stages  of  laminar-turbulent  transition  on  a  flared  cone  in 
Mach  6  flow  have  been  investigated  by  nonlinear  PSE  and  spatial  DNS  techniques.  For  the 
computations,  the  full  effects  of  streamwise  surface  curvature  have  been  taken  into  account. 
These  include  adverse  streamwise  pressure  gradient,  decreased  edge  Mach  number,  a  reversal 
of  boundary-layer  growth,  and  increased  wall  temperature. 

Physical  experiments,  N-factor  studies  based  on  compressible  linear  stability  theory,  non¬ 
linear  PSE  methodology,  and  DNS  all  confirm  that  disturbances  of  very  high  frequency  (on 
the  order  of  230  kHz),  are  responsible  for  transition  from  a  laminar  to  a  turbulent  state  for 
the  flared-cone  configuration  under  study.  Disturbances  in  this  frequency  range  are  classi¬ 
fied  as  “second- mode”  disturbances  in  the  nomenclature  of  Mack  [12].  A  further  indication 
that  the  primary  disturbances  are  of  second-mode  type  is  the  appearance  of  so-called  “rope¬ 
like”  waves  in  schlieren  images  from  the  physical  experiment  or  in  schlieren-like  images  from 
the  numerical  experiment.  These  disturbances  manifest  a  wavelength  that  is  approximately 
twice  the  boundary-layer  displacment  thickness,  in  keeping  with  the  high-frequencey  nature 
of  second-mode  disturbances. 

The  results  obtained  by  nonlinear  PSE  and  DNS  methodologies  show  significant  quan¬ 
titative  disagreement  cis  the  flow  evolves  spatially  from  identical  upstream  initial  states. 
Qualitative  agreement  in  the  growth  rates  and  structures  of  harmonics  persists  quite  far 
downstream.  Despite  these  quantitative  differences,  the  transition  mechanism  identified  as 
viable  by  PSE  methodology  is  confirmed  by  DNS  to  lead  to  a  transitional  state.  The  transi¬ 
tional  state  is  triggered  by  the  interaction  of  a  triad  of  fundamental  disturbance  components 
whose  nonlinear  interactions  generate  all  possible  harmonics  in  (f,m)  wavenumber  space, 
where  integers  I  and  m  refer  to  the  harmonic  with  respect  to  the  fundamental  temporal 
frequency  and  aziumuthal  wavenumber,  respectively.  The  transitional  state  is  characterized 
by  a  decrease  in  the  boundary-layer  shape  factor  relative  to  its  laminar  value  and  increases 
in  both  skin  friction  and  thermal  stress  at  the  wall. 
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The  results  of  the  present  numerical  flared-cone  experiment  differ  dramatically  from  those 
of  Pruett  and  Chang  [16],  who  simulated  laminar  breakdown  on  a  right  circular  cone  without 
a  flared  afterbody.  In  particular,  for  the  present  configuration,  the  dominant  Reynolds  stress 
components,  turbulent  kinetic  energy,  and  turbulent  Mach  number  all  manifest  large  values 
close  to  the  wall,  rather  than  near  the  critical  layer  as  observed  previously.  Moreover,  very 
large  density  fluctuations  occur  both  at  the  wall  and  near  the  critical  layer  for  the  present 
configuration.  For  the  present  problem,  Goertler  modes  of  instability  (i.e.,  harmonics  of  zero 
frequency  of  the  form  (0,m))  are  linearly  unstable.  These  harmonics  manifest  significant 
structure  near  the  wall  and  result  in  earlier  transition  onset  than  is  observed  on  the  straight 
cone,  if  transition  onset  is  defined  on  the  basis  of  skin  friction  or  thermal  stress  at  the  wall. 
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Figure  2:  Boundary-layer  displacement  thickness  versus  surface  arc  length  for  fiared-cone 
model.  Boundary-layer  growth  reverses  at  body  flare. 
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Figure  3:  Normalized  temperature  and  streamwise  velocity  distributions  versus  wall-normal 
coordinate  ^  at  axial  station  x*  =  1.1725  ft. 
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Figure  4:  Normalized  wall-normal  velocity  versus  wall-normal  coordinate  z  at  axial  station 
X*  =  1.1725  ft. 
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Figure  5:  Residuals  of  governing  equations  along  inflow  boundary. 
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respectively,  where  /3*R*  =  15.  ^  7  d  spanwise  wavenuml 
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Figure  7:  Comparison  of  DNS  and  PSE  results  for  selected  harmonics. 
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Figure  8:  Streamwise  evolution  of  selected  mean  quantities. 
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Figure  9:  Streamwise  evolution  of  principal  components  of  Reynolds  stress 
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Figure  10:  Streamwise  evolution  of  turbulent  Mach  number,  turbulent  kinetic  energy,  and 
rms  density  fluctuation. 
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Figure  11:  Amplitudes  of  velocity  fluctuations  of  selected  harmonics  at 
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Figure  13:  Amplitudes  of  temperature  fluctuations  of  selected  harmonics  at  x*  ~  1.05  ft. 
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Figure  14:  Computational  flow  image  of  instantaneous  wall-normal  density  gradient  showing 
“rope-like”  appearance  of  dominant  instability  waves. 
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3.2  DNS  of  Crossflow  Instability  on  a  Supersonic  Swept  Wing 

From  a  computational  point  of  view,  DNS  for  the  elliptic-cone  problem  is  particularly  in¬ 
tensive  because  of  the  absence  of  azimuthal  symmetry.  That  is,  the  computational  domain 
must  include  at  least  one  full  quadrant  of  the  elliptical  cross-section.  In  contrast,  for  the 
circular-cone  geometry  of  Pruett  and  Chang  [16],  symmetry  considerations  in  conjunction 
with  assumptions  about  the  fundamental  azimuthal  wavenumber  =  18)  allowed  re¬ 

striction  of  the  computational  domain  to  just  1/36-th  of  the  circumference  of  the  circular 
cross  section.  Still,  the  simulation  of  Pruett  and  Chang  required  well  in  excess  of  1000  CPU 
hours  on  a  Cray  C90.  Consequently,  to  compute  laminar  breakdown  on  the  elliptic  cone 
could  consume  in  excess  of  10,000  hours  on  existing  supercomputers.  Such  a  computation 
will  be  practical  only  on  the  next  generation  of  supercomputers.  In  the  meantime,  we  believe 
that  it  IS  important  to  develop  the  algorithms  that  will  be  needed  when  the  next  generation 
arrives  and  to  gain  experience  with  the  instability  mechanisms  that  are  likely  to  arise  on 
the  elliptic  cone  and  other  aerospace  configurations.  Moreover,  one  can  continue  to  validate 
other  methodologies  (e.g.,  PSE)  against  DNS  by  considering  linear  and  weakly  nonlinear 
instability  waves,  whose  simulation  is  currently  practicable. 

As  mentioned  previously,  it  is  believed  that  crossflow  instability  will  dominate  on  the 
elliptic  cone,  but  no  experience  exists  in  simulating  crossflow  instability  for  high-speed  flows 
other  than  the  PSE  calculations  of  Chang  et  al.  [5].  By  considering  a  different  configuration 
than  the  elliptic  cone,  namely  the  infinitely  long  swept  wing,  we  have  been  able  to  gain 
experience  with  the  simulation  of  crossflow  instability  in  a  supersonic  boundary-layer  flow. 
The  configuration  considered  is  the  77-degree  swept  wing  of  Cattafesta  et  al.  [3],  which  is 
undergoing  stability  and  transition  tests  in  NASA  Langley’s  Mach  3.5  quite  wind  tunnel. 
Near  the  aft  station  of  the  wing,  cross  sections  axe  similar  and  pressure  isobars  are  nearly 
parallel  to  the  leading  edge.  Consequently,  one  can  approximate  the  flow  as  quasi-three- 
dimensional;  that  is,  the  base  flow  varies  only  in  the  streamwise  and  wall-normal  directions 
and  not  in  the  spanwise  direction  (parallel  to  the  leading  edge).  This  approximation  is 
equivalent  to  assuming  that  the  wing  is  infinitely  long  in  span  and  permits  a  spanwise 
periodicity  assumption  for  the  perturbed  flow.  As  a  result  of  the  particular  periodicity 
assumption,  DNS  can  be  performed  with  reasonable  CPU  requirements. 

This  work  was  begun  under  NASA  contract  NASl-19656  and  was  brought  to  the  current 
state  under  the  present  AFOSR  grant.  The  production  run  was  accomplished  under  the 
auspices  of  AFOSR  on  the  Waterways  Experiment  Station  (WES)  Cray  C90  at  Vicksburg, 
MS.  Progress  to  date  is  documented  thoroughly  in  Attachment  2,  a  NASA  Contractor  Report 
entitled  “Simulation  of  Crossflow  Instability  on  a  Supersonic  Highly  Swept  Wing,”  and  the 
reader  is  referred  to  the  attachment  for  further  details.  At  present  we  have  considered  only 
stationary  crossflow  modes;  traveling  crossflow  modes  also  need  to  be  considered. 
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3.3  Setup  of  Elliptic-Cone  Problem 


To  avoid  confusion  with  the  flared-cone  work  presented  earlier,  please  note  that  different 
coordinate  notation  is  used  for  the  elliptic-cone  problem.  Here  x,  t/,  and  2:  denote  major 
axis,  minor  axis,  and  axial  coordinates  respectively, 

•  The  few  available  documents  (Refs.  [7],  [11],  and  [19])  on  the  elliptic-cone  problem 
were  reviewed.  It  was  realized  that,  in  essence,  the  elliptic-cone  configuration  can  be 
considered  a  “fat”  and  symmetrical  highly  swept  wing,  and  if  the  eccentricity  ratio  is 
moderately  large,  the  flow  will  be  dominated  by  crossflow  instability.  This  realization 
precipitated  the  swept-wing  study  mentioned  in  the  previous  subsection. 

•  Some  thought  was  given  to  determination  of  the  laminar  base  state,  which  may  turn 
out  to  be  the  most  difficult  task  for  the  elliptic-cone  problem.  Unlike  all  other  prob¬ 
lems  addressed  previously  by  this  P.I.,  the  elliptic-cone  problem  presents  ultimately 
the  difficulty  of  a  fully  three-dimensional  base  flow.  Two  possible  options  currently 
exist  to  determine  the  basic  state.  Both  Lyttle  and  Reed  [11]  and  Huang  et  al.  [7] 
exploit  parabolized  Navier-Stokes  (PNS)  technology  to  determine  the  base  flow.  How¬ 
ever,  Lyttle  and  Reed  state:  “The  treatment  of  the  subsonic  pressure  gradient  terms 
introduces  errors  that  make  the  calculated  basic  state  unsuitable  for  stability  analysis.” 
In  our  experience,  PNS  calculations  often  require  numerical  dissipation  that  may  ax- 
tificially  but  significantly  alter  the  stability  of  the  base  state.  Moreover,  well-resolved 
PNS  computations  can  be  expensive  computationally,  requiring  on  the  order  of  100 
CPU  hours  on  a  supercomputer.  Consequently,  we  prefer  at  the  present  time  to  use 
an  Euler/boundary-layer  approach  whereby  the  outer  flow  is  determined  by  solution  of 
the  Euler  equations  and  the  inner  solution  is  obtained  from  the  boundary-layer  equa¬ 
tions.  Fortunately,  for  conical  flows,  including  flow  over  an  elliptic  cone,  the  inviscid 
flow  does  not  vary  along  cone  generators,  as  will  be  verified  later.  Once  again,  this 
permits  a  simplifying  assumption  by  which  the  flow  can  be  approximated  as  quasi- 
three-dimensional.  It  is  not  altogether  certain  that  this  approach  will  work,  however, 
because  the  conical-flow  boundary- layer  solution  will  not  automatically  satisfy  symme¬ 
try  conditions.  Our  intuition  is  that  the  approach  will  work  for  the  case  of  zero  angle 
of  attack  but  that  viscous-inviscid  interaction  will  become  too  strong  when  the  cone  is 
at  non-zero  incidence  angle  to  the  flow. 

•  An  elliptic-cone  configuration  was  defined  based  on  the  model  used  by  Schneider  [19]  for 

his  Ludwieg-tube  experiments;  nominal  freestream  conditions  were  also  obtained  from 
Schneider’s  experiment.  The  model  is  5  inches  long  with  a  sharp  apex  and  elliptical 
cross  sections.  The  eccentricity  ratio  is  4:1,  and  the  major-axis  semi-vertex  angle  is  17.5 
degrees.  At  a  freestream  unit  Reynolds  number  of  Re\  =  1,431,070,  the  freestream 
conditions  are  M^o  =  4.0,  =  13.66  psf.,  and  =  128.7  Rankine,  for  a  ratio  of 

specific  heats,  7  =  1.4. 
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•  The  Euler  equations  were  solved  for  the  configuration  and  flow  conditions  above. 
Specifically,  the  marching  Euler  code  “ARROW”  of  Salas  [18]  was  used  to  obtain  the 
inviscid  solution.  The  pressure  coefficient  so  obtained  is  shown  versus  the  scaled  minor- 
axis  coordinate  y  in  Fig.  15,  for  several  streamwise  stations  2.  Note  that,  in  scaled 
coordinates,  all  profiles  collapse  onto  one  as  the  solution  is  marched  downstream;  thus, 
the  conical-flow  assumption  is  verified  for  this  configuration.  Since  the  starting  condi¬ 
tion  for  ARROW  assumes  the  cross-section  is  circular,  the  initial  profiles  are  in  error 
but  do  not  significantly  affect  the  2-asymptotic  state. 

•  The  governing  boundary-layer  equations  for  conical  flow  were  derived  so  that  the  spec¬ 
trally  accurate  boundary-layer  code  of  Pruett  and  Streett  [13]  could  be  modified  to 
incorporate  conical  flows.  The  modification  to  the  code  has  not  yet  been  accomplished. 

•  An  interpretation  routine  was  written  to  convert  the  output  of  the  inviscid  solution  into 
edge  conditions  for  the  boundary-layer  code.  The  boundary-layer  equations  exploit  a 
body-fitted  coordinate  system  for  which  coordinates  x\  y',  and  z'  specify  the  distance 
along  the  cone  generators,  the  arc  length  perpendicular  to  the  generator,  and  the  wall- 
normal  distance,  respectively.  The  corresponding  edge  velocities,  u*,  v*,  and  w*,  are 
shown  in  Figs.  16,  17,  and  18,  respectively.  By  definition,  the  normal  velocity  lo* 
should  vanish  for  inviscid  flow;  hence,  the  small  values  shown  in  Fig.  18  serve  purely 
as  a  check  that  the  inviscid  solution  hcis  been  properly  computed  and/or  interpreted. 

•  Some  thought  was  given  to  derivation  of  the  compressible  Navier-Stokes  equations  in  a 
body-fitted  coordinate  system  for  the  elliptic-cone  configuration,  as  was  done  by  Pruett 
et  al.  [15]  for  the  cone  and/or  flared-cone  configurations.  Unfortunately,  the  absence 
of  azimuthal  symmetry  adds  a  quantum  leap  in  complexity  to  the  metric  terms  for  the 
elliptic-cone  problem.  Thus,  it  is  our  judgement  that  the  elliptic-cone  problem  is  best 
treated  by  resorting  to  generalized  curvilinear  coordinates.  From  a  computational  point 
of  view,  generalized  coordinates  will  likely  increase  CPU  and  memory  requirements  each 
by  a  factor  of  two  relative  to  the  current  baseline  code.  The  advantage,  however,  is 
that  configurations  other  than  the  elliptic  cone  could  be  treated  readily  in  a  generalized 
coordinate  system. 

•  Additional  thought  was  given  to  optimal  parameterization  of  the  azimuthal  coordinate 
for  the  elliptic-cone  configuration.  Huang  et  al.  [7]  parameterize  as  if  the  cross-section 
were  circular;  that  is,  based  on  the  angle  formed  by  rays  from  the  center  to  points  on 
the  surface.  This  parameterization  has  the  unfortunate  trait  of  placing  points  further 
apart  in  regions  where  curvature  is  largest.  We  favor  parameterization  based  on  the 
“clock  angle”  ^  formed  by  the  surface  normal.  Fig.  19  contrasts  the  two  parameteriza- 
tions  for  the  edge  Mach  number  Mg.  Note  that  the  Mach  number  distribution  is  much 
smoother  in  the  preferred  parameterization. 
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Figure  15:  Inviscid  surface  pressure  coefficient  versus  scaled  minor-axis  coordinate  at  three 
axial  stations  obtained  with  ARROW  code  of  Salas.  Collapse  of  all  three  profiles  onto  one 
curve  indicates  that  flow  is  conical  on  elliptic  cone.  (Length  of  major  semi-axis  is  Cn.) 
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17.  Contours  of  constant  disturbance  density  in  computational  plane  nor¬ 
mal  to  both  leading  edge  and  wing  surface  showing  acoustic  energy 
generated  by  strong  nonlinearity  in  wave-packet  region. 


parameterization  angle 


Figure  19:  Comparison  of  two  possible  parameterizations  of  azimuthal  coordinate  for  elliptic- 
cone  problem.  Proposed  clock-angle  parameterization  results  in  better  distribution  of  grid 
points. 
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4  Personnel 


This  grant  has  supported  the  research  efforts  of  the  Principal  Investigator,  Dr.  C.  David 
Pruett,  Dept,  of  Applied  Science,  The  College  of  William  and  Mary,  Williamsburg,  VA. 
In  addition.  Dr.  Chau-Lyan  Chang  of  High  Technology  Corporation,  Hampton,  VA,  has 
contributed  substantially  to  this  effort  without  remuneration  from  AFOSR. 


5  Publications 

•  C.  D.  Pruett  and  C.-L.  Chang,  “Spatial  Direct  Numerical  Simulation  of  High-Speed 
Boundary-Layers  Flows-Part  II:  Transition  on  a  Cone  in  Mach  8  Flow,”  Theoretical 
and  Computational  Fluid  Dynamics,  Vol.  7,  No.  5,  1995,  pp.  397-424. 

•  C.  D.  Pruett  and  T.  A.  Zang,  “On  Simulation  and  Analysis  of  Instability  and  Transition 
in  High-Speed  Boundary-Layer  Flows,”  Computing  Systems  in  Engineering,  to  appear. 

•  C.  D.  Pruett,  “Spatial  Direct  Numerical  Simulation  of  Transitioning  High-Speed  Flows” 
(Invited  Paper),  presented  to  the  Second  Symposium  on  Transitioning  and  Turbu¬ 
lent  Compressible  Flows  at  the  ASME/JSME  Fluids  Engineering  Annual  Conference, 
Hilton  Head,  SC,  13-18  August,  1995. 


6  Interactions  /  Transitions 

•  Addressed  the  Virginia  Consortium  of  Engineering  and  Science  Universities  (VCES) 
Seminar  in  Hampton,  VA,  on  10  March,  1995,  on  the  topic  of  “Direct  Numerical 
Simulation  of  Transitioning  High-Speed  Flows.” 

•  Visited  AFOSR  at  Bolling  AFB,  DC,  on  24  July,  1995,  to  discuss  current  and  future 
directions  of  this  effort  with  the  technical  monitor.  Dr.  Leonidas  Sakell. 

•  Addressed  the  Second  Symposium  on  Transitional  and  Turbulent  Compressible  Flows 
at  the  ASME/JSME  Fluids  Engineering  Annual  Conference,  Hilton  Head,  SC,  13-18 
August,  1995.  (See  “Publications.”) 


7  Discoveries,  Inventions,  or  Patents 

A  qualitative  difference  exists  between  transitioning  hypersonic  boundary-layer  flows  on  ax- 
isymmetric  straight  and  flared  cones.  For  the  case  of  the  straight  cone,  Reynolds  stresses  are 
largest  near  the  critical  layer,  spreading  toward  the  wall  only  through  nonlinear  effects.  For 
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the  flared  cone,  Goertler  modes,  associated  with  concave  streamwise  curvature,  are  linearly 
unstable  and  introduce  large  fluctuations  close  to  the  wall  in  the  the  relatively  early  stages  of 
transition.  As  a  result,  Reynolds  stresses  peak  close  to  the  wall.  From  an  engineering-design 
point  of  view,  transition  onset,  as  determined  by  increased  wall  shear  and  thermal  transfer, 
will  be  observed  sooner  in  the  flared-cone  case.  However,  the  transition  process  is  relatively 
gradual  and  the  transition  zone  may  be  relatively  long  for  the  flared  cone. 


8  Honors  and  Awards 


Invited  to  address  the  Symposium  on  Transitional  and  Turbulent  Compressible  Flows  at  the 
ASME/JSME  Fluids  Engineering  Annual  Conference,  Hilton  Head,  SC,  13-18  August,  1995. 
(See  “Interactions”  and  “Publications.”) 


9  Attachments 

1.  Paper  by  Pruett,  Zang,  Chang,  and  Carpenter  entitled  “Spatial  Direct  Numerical  Sim¬ 
ulation  of  High-Speed  Boundary-Layer  Flows-Part  I:  Algorithmic  Considerations  and 
Validation.” 

2.  NASA  Contractor  Report  by  Pruett  entitled  “Simulation  of  Crossflow  Instability  on  a 
Supersonic  Highly  Swept  Wing.” 
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Abstract.  A  highly  accurate  algorithm  for  the  direct  numerical  simulation  (DNS)  of  spatially 
evolving  high-speed  boundary-layer  flows  is  described  in  detail  and  is  carefully  validated.  To 
represent  the  evolution  of  instability  waves  faithfully,  the  fully  explicit  scheme  relies  on  non- 
dissipative  high-order  compact-difference  and  spectral  collocation  methods.  Several  physical, 
mathematical,  and  practical  issues  relevant  to  the  simulation  of  high-speed  transitional  flows  are 
discussed.  In  particular,  careful  attention  is  paid  to  the  implementation  of  inflow,  outflow,  and 
far-field  boundary  conditions.  Four  validation  cases  are  presented,  in  which  comparisons  are 
made  between  DNS  results  and  results  obtained  from  either  compressible  linear  stability  theory  or 
from  the  parabolized  stability  equation  (PSE)  method,  the  latter  of  which  is  valid  for  nonparallel 
flows  and  moderately  nonlinear  disturbance  amplitudes.  The  first  three  test  cases  consider  the 
propagation  of  two-dimensional  second-mode  disturbances  in  Mach  4.5  flat-plate  boundary-layer 
flows.  The  final  test  case  considers  the  evolution  of  a  pair  of  oblique  second-mode  disturbances  in 
a  Mach  6,8  flow  along  a  sharp  cone.  The  agreement  between  the  fundamentally  different  PSE  and 
DNS  approaches  is  remarkable  for  the  test  cases  presented. 


L  Introduction 

A  worthy  “grand  challenge”  for  the  computational  boundary-layer-transition  community  is  the  accu¬ 
rate  direct  numerical  simulation  (DNS)  of  the  complete  laminar-turbulent  transition  process  in  a 
spatially  evolving  high-speed  boundary-layer  flow.  Even  for  such  simple  geometries  as  the  flat  plate  or 
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sharp  cone,  this  remains  a  daunting  goal.  It  was  only  in  the  late  1980s  that  the  corresponding 
challenge  for  temporally  evolving  incompressible  flow  was  met  by  Gilbert  and  Kleiser  (1990)  for  the 
relatively  simple  problem  of  channel-flow  transition.  For  spatially  evolving  compressible  boundary- 
layer  flows,  the  recent  landmark  simulations  of  Rai  and  Moin  (1991)  and  Thumm  (1991)  have  come 
closest  to  the  realization  of  this  goal.  In  the  former  simulation  a  fully  turbulent  state  was  attained 
from  an  initially  laminar  state  subject  to  high-amplitude  random  forcing  imposed  in  the  free  stream. 
However,  although  the  algorithm  was  designed  for  compressible  flow,  the  low  subsonic  Mach  number 
(Mach  0.1)  of  the  numerical  experiment  guaranteed  that  the  flow  was  essentially  incompressible  in 
behavior.  The  latter  computation  of  Thumm  (1991)  considered  a  relatively  low-speed  supersonic  flow 
(Mach  1.6),  and,  while  simulating  highly  nonlinear  stages  of  transition,  it  did  not  proceed  into  a 
fully  turbulent  regime.  (For  a  thorough  background  of  DNS  for  transitional  incompressible  and 
compressible  wall-bounded  flows,  including  a  discussion  of  the  temporal  and  spatial  problems,  see 
Kleiser  and  Zang  (1991).) 

The  first  tentative  steps  were  taken  in  the  use  of  DNS  to  investigate  transition  to  turbulence  in 
supersonic,  wall-bounded  flows  in  the  mid  to  late  1980s.  Bayliss  et  al.  (1985)  presented  the  first  DNS 
results  for  supersonic  boundary-layer  flow  along  a  flat  plate.  These  results  were  for  spatially  evolving, 
but  two-dimensional,  flow.  The  first  three-dimensional  DNS  of  a  perturbed  high-speed  (Mach  4.5) 
flat-plate  boundary-layer  flow  was  accomplished  by  Erlebacher  and  Hussaini  (1990),  This  numerical 
experiment  used  temporal  DNS  to  examine  boundary-layer  stability,  but  stopped  far  short  of 
attaining  a  transitional  state. 

Recently,  due  partly  to  increased  supercomputer  capacity,  there  have  been  many  noteworthy 
three-dimensional  simulations  of  compressible  wall-bounded  flows.  Among  these  are  temporal  simula¬ 
tions  by  Normand  and  Lesieur  (1992),  Pruett  and  Zang  (1992),  Dinavahi  and  Pruett  (1993),  and 
Adams  and  Kleiser  (1993);  and  spatial  simulations  by  Thumm  et  al  (1990),  Maestrello  et  al  (1991), 
Thumm  (1991),  Normand  and  Lesieur  (1992),  Eissler  and  Bestek  (1993),  Ng  and  Zang  (1993),  and 
Pruett  and  Chang  (1993).  Among  these,  the  temporal  simulation  of  Dinavahi  and  Pruett  (1993)  is 
unique  in  attaining  a  well-resolved  fully  turbulent  state  without  recourse  to  modeling. 

To  date,  with  the  exception  of  Rai  and  Moin  (1991),  virtually  all  of  the  numerical  experiments  cited 
have  simulated  controlled  rather  than  natural  instability  processes.  In  a  controlled  experiment, 
instability  waves  of  a  particular  wavelength  (temporal)  or  frequency  (spatial)  are  excited  by  imposed 
forcing.  In  contrast,  in  natural  transition  the  input  is  random,  and  the  flow  itself  selects  the  preferred 
instability  modes.  A  few  of  the  cited  simulations  are  hybrid  in  the  sense  that  the  primary  instability 
wave  is  imposed,  whereas  secondary  instability  is  triggered  by  low-level  noise.  A  distinguishing  feature 
of  high-speed  boundary-layer  flows  is  that  multiple  primary  instability  modes  can  coexist  (Mack, 
1984).  The  viscous  first-mode  instability,  the  counterpart  of  the  Tollmien-Schlichting  (TS)  wave  in 
incompressible  flow,  predominates  in  low-speed  compressible  flows.  In  hypersonic  boundary-layer 
flows  second-mode  instabilities  arise,  which  are  acoustic  in  nature,  and  which  eventually  predominate 
as  the  Mach  number  increases.  Thus  far,  all  simulations  of  instability  and  transition  in  compressible 
flows  have  focused  on  first-  and  second-mode  instabilities  and  their  associated  secondary  instabilities, 
rather  than  on  crossflow  or  Goertler  modes. 

Several  observations  gleaned  from  temporal  DNS  are  worthy  of  note.  Normand  and  Lesieur  (1992) 
performed  both  DNS  of  a  low-speed  (Mach  0.5)  compressible  flat-plate  flow  and  large-eddy  simula¬ 
tion  of  a  high-speed  (Mach  5)  flow.  They  observed  transition  to  occur  by  means  of  fundamental 
secondary  instability  in  the  former  case  and  subharmonic  secondary  instability  in  the  latter  case.  Their 
findings  are  consistent  with  those  of  Pruett  and  Zang  (1992),  Dinavahi  and  Pruett  (1993),  and  Adams 
and  Kleiser  (1993).  Pruett  and  Zang  (1992)  and  Dinavahi  and  Pruett  (1993)  considered  the  case  or 
Mach  4.5  flow  along  a  hollow  cylinder  (the  axisymmetric  analog  of  a  flat-plate  boundary  layer)  and 
observed  subharmonic  secondary  instability  triggered  by  second-mode  primary  instability  to  lead  to 
transition.  Adams  and  Kleiser  (1993)  performed  a  similar  computation  for  a  Mach  4.5  fat  plate  using 
random  noise  to  trigger  secondary  instability  from  a  base  state  perturbed  by  a  second-mode  distur¬ 
bance.  The  secondary  instabilities  they  observe  are  of  subharmonic  type  and  agree  in  growth  rate  and 
structure  with  the  predictions  of  temporal  secondary  instability  theory  (Ng  and  Erlebacher,  1992). 
Moreover,  despite  the  difference  in  geometry,  their  results  are  in  qualitative  agreement  with  those  of 
Pruett  and  Zang  (1992).  A  potential  weakness  of  each  of  these  simulations,  however,  is  the  failure  to 
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Figure  1.  Amplitudes  of  components  of  oblique  second-mode  Q.O  0.5  1.0  1.5  2.0 
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account  for  growth  of  the  boundary  layer,  a  somewhat  ambiguous  task  in  the  context  of  temporal 
theory.  Recent  experimental  results  (Stetson  and  Kimmel,  1993)  and  numerical  results  obtained 
from  the  parabolized-stability-equation  (PSE)  method  (Chang,  1993)  suggest  that  subharmonic  sec¬ 
ondary  instability  may  not  be  the  preferred  path  to  transition  in  a  growing  high-speed  boundary 
layer.  The  issue  can  probably  only  be  resolved  by  spatial  DNS,  which  incorporates  the  evolution  of 
the  boundary  layer. 

The  temporal  simulations  cited  above  required  upward  of  10®  grid  points  and  consumed  hundreds 
of  Cray  Y-MP  CPU  hours.  There  are  several  reasons  for  the  great  expense  of  computations  of 
high-speed  compressible  flow  relative  to  simulations  of  incompressible  or  subsonic  flows: 

1.  The  time  discretizations  in  the  compressible  cases  were  fully  explicit.  In  many  instances  the 
allowable  time  step  was  limited  by  the  viscous  stability  condition  rather  than  by  the  advection 
condition.  For  incompressible  simulations  the  viscous  stability  limit  is  usually  absent  due  to  the 
conventional  implicit  treatment  of  the  viscous  terms. 

2.  The  second-mode  disturbances  associated  with  high-speed  transitional  flows  have  a  double- 
peaked  structure  (Figure  1)  with  amplitude  peaks  occurring  both  near  the  wall  and  the  critical 
layer  (z  a  1  in  Figure  1).  In  contrast  to  low-speed  flows,  at  high  speeds  the  critical  layer  lies  far 
(approximately  one  displacement  thickness)  from  the  wall,  necessitating  concentrations  of  grid 
points  in  both  regions  of  strong  gradients. 

3.  At  high  speeds  the  growth  rates  of  both  the  primary  and  secondary  instabilities  are  much  slower 
than  for  low-speed  flows.  This  requires  much  longer  time  integrations. 

4.  In  contrast  to  DNS  of  incompressible  flow,  flow-field  oscillations  due  to  inadequate  resolution 
are  potentially  fatal  in  the  compressible  case  since  spurious  negative  densities,  pressures,  and/or 
temperatures  can  arise. 

Relative  to  temporal  DNS,  spatial  DNS  is  yet  more  computationally  demanding,  primarily  because 
of  the  greater  length  of  the  computational  domain.  Nonetheless,  the  computational  boundary-layer 
transition  community  has  recently  produced  several  spatial  simulations  of  particular  importance.  The 
simulations  of  Thumm  et  al.  (1990)  examine  the  very  early  stages  of  secondary  instabilities  of  both 
fundamental  and  subharmonic  type  in  Mach  1.6  boundary-layer  flows  along  a  fiat  plate.  Subse¬ 
quently,  Thumm  (1991)  turned  his  attention  to  investigating  a  new  “oblique-mode”  transition  mecha¬ 
nism,  again  for  Mach  1.6  fiat-plate  flow.  Conventionally,  transition  is  triggered  by  secondary  instabil¬ 
ity  originating  from  a  finite-amplitude  two-dimensional  primary  disturbance  and  a  pair  (or  array)  of 
low-amplitude  oblique  secondary  disturbances.  For  Mach  1.6  boundary-layer  flow,  however,  the  most 
amplified  linear  instability  turns  out  to  be  oblique  and  not  two-dimensional,  in  contrast  to  the  case 
for  either  incompressible  flow  or  flow  at  very  high  Mach  numbers.  This  suggests  the  possibility  of 
triggering  transition  simply  by  the  nonlinear  interactions  of  a  pair  of  symmetric  unstable  primary 
modes.  Thumm  (1991)  presented  an  extensive  set  of  results,  comparing  fundamental-,  subharmonic-. 
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and  oblique-breakdown  scenarios.  The  simulations,  which  were  carried  into  the  early  nonlinear  stages 
of  laminar  breakdown,  have  recently  been  summarized  by  Fasel  et  al  (1993).  Maestrello  et  al 
(1991)  performed  three-dimensional  spatial  simulations  of  Mach  4.5  flows  excited  simultaneously  by  a 
two-dimensional  second-mode  disturbance  and  a  single  oblique  disturbance  of  the  same  frequency. 
They  computed  to  the  early  nonlinear  stages  of  transition  and  demonstrated  significant  interactions 
between  the  forced  modes,  which  then  generated  other  rapidly  growing  instability  modes.  Eissler  and 
Bestek  (1993)  performed  several  simulations  at  Mach  4.8  using  periodic  suction  and  blowing  at  the 
wall  to  excite  instability  waves.  In  addition  to  the  excitation  of  the  expected  oblique  second  mode, 
their  results  showed  excitation  of  an  additional  “viscous”  mode  of  the  same  frequency  but  of  dilferent 
wavelength.  The  evolution  of  these  modes  was  tracked  from  the  first-mode  region  into  the  second¬ 
mode  region. 

Of  these  simulations,  the  typical  spatial  accuracy  was  fourth-order,  occasionally  with  a  spectrally 
accurate  Fourier  collocation  method  in  the  (periodic)  spanwise  direction.  In  most  of  these  spatial 
simulations  the  spanwise  resolution  was  modest,  in  some  cases  admitting  just  a  single  oblique  mode. 
Although  Thumm  (1991)  used  up  to  17  grid  points  in  the  spanwise  direction,  and  Maestrello 
et  al.  (1991)  used  up  to  32  spanwise  grid  points,  such  resolution  is  believed  to  be  insufficient  for  the 
later  stages  of  transition,  in  which  there  is  typically  an  explosive  broadening  in  wave  space  of  the 
energy  spectrum  in  the  spanwise  direction.  The  simulations  at  the  higher  Mach  numbers  call  attention 
to  the  difficulties  (and  impracticality)  of  simulating  the  entire  laminar-turbulent  transition  process 
when  the  dominant  second-mode  instabilities  are  of  high  frequency  but  slow  growth.  To  begin  from  a 
linearly  perturbed  laminar  state,  the  computational  domain  must  be  extremely  long  in  streamwise 
extent  (relative  to  the  disturbance  wavelength).  Alternately,  the  forcing  amplitudes  must  be  quite 
large,  in  which  case  nonlinear  interactions  are  significant  at  the  inflow  boundary  and  consistent  inflow 
conditions  are  difficult  to  obtain. 

This  work  focuses  on  the  development  of  an  algorithm  and  additional  procedures  for  which  spatial 
DNS  of  such  challenging  high-speed  transition  problems  is  feasible.  With  regard  to  the  algorithm,  its 
origins  can  be  found  in  Ng  and  Zang  (1993),  who  compared  spatial  DNS  with  secondary  instability 
predictions  for  Mach  1.6  and  Mach  6.8  flows  in  artificial,  quasi-parallel  boundary-layer  flows. 
Their  work  validated  both  the  basic  spatial  DNS  code  and  the  spatial  secondary  instability  code,  but 
it  was  not  an  investigation  into  the  physics  of  transition.  Their  work,  however,  formed  the  starting 
point  for  the  algorithm  developed  by  Pruett  and  Chang  (1993)  for  the  spatial  DNS  computations  and 
the  results  reported  in  this  paper.  Pruett  and  Chang  (1993)  made  detailed  comparisons  of  spatial  DNS 
and  PSE  for  two-dimensional  linear  and  nonlinear  second-mode  disturbances  in  Mach  4,5  flat-plate 
boundary  layers.  They  achieved  remarkable  agreement  between  the  methods  for  these  cases;  however, 
they  found  that  the  numerical  methods  utilized  in  the  DNS  had  to  be  extremely  refined  in  order  to 
obtain  accurate  results.  Careful  attention  was  given  to  many  details  often  overlooked.  Among  the 
issues  addressed  were:  obtaining  self-consistent  inflow  conditions  for  nonparallel  flows,  the  necessity  of 
filtering  to  suppress  spurious  high-frequency  modes  and  boundary  reflections,  obtaining  reliable 
estimates  for  the  numerical  stability  limits  on  the  time  step  as  a  function  of  the  precise  spatial 
discretization,  and  the  imposition  of  nonreflecting  far-field  and  outflow  boundary  conditions. 

Our  present  position  is  that  the  “grand  challenge”  of  simulating  transition  in  high-speed  boundary- 
layer  flows  is  best  met  by  a  combination  of  numerical  tools.  The  PSE  method  is  much  more  efficient 
than  DNS  for  computing  the  early  linear  and  weakly  nonlinear  stages  of  the  transition  process.  Thus, 
we  advocate  the  use  of  PSE  to  furnish  accurate  and  consistent  inflow  conditions  to  a  DNS  that 
commences  near  the  highly  nonlinear  laminar-breakdown  stage.  The  primary  objective  of  this  paper  is 
to  provide  thorough  validation  and  documentation  of  the  DNS  algorithm.  In  the  cross-validation  of 
DNS  with  PSE,  however,  we  demonstrate  convincingly  the  fidelity  of  the  PSE  method  within  the 
appropriate  flow  regimes. 

An  alternative,  but  still  developing  approach,  is  that  of  large-eddy  simulation  (LES).  In  this  method 
large-scale  structures  are  resolved,  but  a  model  is  employed  for  the  subgrid-scale  fluctuations.  The 
LES  approach  is  fairly  well  developed  for  incompressible  flows  (see  Kleiser  and  Zang,  1991),  but  only 
a  few  results  are  available  for  supersonic  flows  (see  Normand  and  Lesieur,  1992;  Zang  et  al, 
1992;  El-Hady  et  al,  1993).  Because  the  cost  of  PSE  increases  much  faster  than  that  of  DNS  as  the 
number  of  spanwise  modes  is  increased,  PSE  is  limited  in  practice  to  investigations  of  narrow-band 
forcing.  At  present,  LES  appears  to  be  the  only  hope  for  complete  simulations  of  natural  (broad-band) 
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transition.  In  the  same  way  that  the  spatial  DNS  described  herein  has  served  to  validate  PSE,  it  is 
hoped  that  it  will  also  be  useful  in  validating  LES  approaches. 

This  paper  consists  of  a  reprise  of  the  Pruett  and  Chang  (1993)  results,  with  an  emphasis  on  the 
subtle  numerical  issues  that  are  essential  for  accurate  spatial  DNS  of  high-speed  boundary-layer  flows. 
Their  previous  work  is  also  extended  to  the  simulation  of  three-dimensional  second-mode  distur¬ 
bances  propagating  in  the  boundary  layer  of  a  sharp  cone  immersed  in  a  Mach  8  flow.  The 
motivation  for  this  example  is  the  stability  experiment  by  Stetson  et  al  (1983),  in  which  the 
second-mode  disturbance  was  observed  as  the  dominant  instability  in  such  a  flow.  The  ultimate 
objective  of  this  work,  and  the  subject  of  the  sequel  to  this  paper,  is  the  combined  PSE/DNS 
investigation  of  the  complete  laminar-turbulent  transition  process  in  a  hypersonic  boundary-layer  flow 
along  a  cone. 

In  the  next  section  the  three-dimensional  compressible  Navier- Stokes  equations  are  presented  for 
flow  along  a  two-dimensional  or  an  axisymmetric  body.  The  fully  explicit  numerical  method,  which 
combines  high-order  compact-difference  and  spectral  collocation  methods,  is  presented  in  detail  in 
Section  3.  Algorithm  details  are  discussed  in  Section  4.  In  Section  5  compressible  linear  stability 
theory  (LST)  and  the  PSE  method  are  discussed  very  briefly,  and  the  DNS  code  is  thoroughly 
validated  against  these  yardsticks.  Three  of  the  four  validation  cases  examine  the  evolution  of 
two-dimensional  second-mode  instability  waves  in  a  Mach  4.5  flat-plate  boundary-layer  flow.  Both 
linear  and  nonlinear  disturbance  amplitudes  are  considered.  The  fourth  case  examines  the  evolution  of 
a  pair  of  oblique  second-mode  disturbances  on  a  sharp  T  half-angle  cone  in  a  Mach  6.8  (free-stream 
Mach  8)  flow.  Conclusions  which  relate  both  to  the  numerical  method  and  to  transition  physics  are 
presented  in  the  final  section. 


2.  Governing  Equations 


Consider  the  body-fitted  orthogonal  coordinate  system  x  =  [x,  0,  z]^  on  an  axisymmetric  body  as 
shown  in  Figure  2,  where  x  is  the  arc  length  along  the  body,  0  is  the  azimuthal  angle,  and  z  is  the 
coordinate  normal  to  the  body.  Associated  with  x  is  the  fundamental  metric  tensor  which  has  the 
nonzero  components 

0U=S^  g22  =  r^,  033  =  U  (1) 

where 


r(x,  z)  =  R  z  cos  (p, 

s{x,  z)  =  1  -  z^, 
ax 


(2) 


and  where  R{x)  and  ^(x)  are  the  body  radius  and  the  angle  of  the  surface  tangent  in  the  plane  of 
symmetry,  respectively.  For  convenience,  we  define  the  following  partial  differential  operators,  which 
incorporate  the  metric  quantities  r  and  s: 


Dy  =  -—, 

s  ox 


0„  — 


DSu 


1  du 

rW 


du 


Diu  = 


1  d{ru) 
rs  dx 


,  _  1  5(rsu) 


D}u  = 


rs  dz 


(3) 


Figure  2.  Body-fitted  coordinate  system  on  axisymmetric 
body. 


54 


C.D.  Pruett,  T.A.  Zang,  Chau-Lyan  Chang,  and  M.H.  Carpenter 


In  the  coordinate  system  of  Figure  2  and  in  the  terminology  of  (3),  the  dimensionless  compressible 
Navier-Stokes  equations  assume  the  following  conservative  form: 

^  +  Z)iE  +  D9«F  +  DiG  +  H  =  0,  (4) 

ct 

where  Q  =  [p,  pu,  pv,  pw,  E.f  is  the  fluid  state  vector  of  conserved  quantities;  p  is  the  density; 
u,  V,  and  w  are  velocities  in  the  e^,  e^,  and  directions,  respectively;  and  £,  is  the  total  energy.  For 
later  use,  we  also  define  the  vector  of  primitive  variables  U  =  [p,  u,  v,  w,  pY,  and,  for  state  and  flux 
vectors  in  general,  Q  =  [_Qo,  Qu  Qi,  Qs,  Q^li\  etc.  Vectors  E,  F,  G,  and  H  are  defined  as  follows: 

pu 

pm  —  Til  ”  P 

E  =  puv  —  Ti2  , 

puw  -  Ti3 

(£,  +  p)u  -  UTii  -  VX^2  - 


pVU  -  T21 
pvv  —  T22  P 

puw  -  T23 

+  p)v  -  MT21  -  VT22  -  WT23  -  ^2 
pw 

pwu  -  T31 
pwv  -  T32 
pww  -  T33  -  p 


+  p)w  -  MT31  -  UT32  -  WT3 


-F2  -  (i)»£3 


sin(p  coscp 

- 

r  r 


where  p  is  thermodynamic  pressure. 


-F2  +  (Z)»£i 


{y  -  l)M^Re^'‘ 
^  (7  - 


u  =  _ po  T 

(7  -  l)M^Re  " 

are  the  heat  flux  components,  T  is  the  temperature,  k  is  the  thermal  conductivity, 

is  the  stress  tensor,  p  is  the  dynamic  viscosity,  is  the  Kronecker  delta,  d  =  D^u  +  DgV  +  D^w  is  the 
divergence  of  the  velocity  (dilatation),  and  Uy  is  the  symmetric  rate-of-deformation  tensor  with 
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components 


=  D^u  — 


^  usm  (D  W  cos  (D 

_  no*,  j _ 1.  j - — 


^22  ~  ^0  ^ 


(J32  _  w, 

^12  “  ^21  “  2(^6 y 


V  sincp 


<^13  =  <^31  =  2(^>  +  -  ^^X<P\ 

.  f  n  n  ^  cos  q>\ 

(723  =  ^32  =  +  - - - j- 

With  the  exception  of  (1),  which  refers  to  a  co variant  tensor,  all  other  equations  refer  to  physical 
vector  or  tensor  components.  Thus,  for  convenience  and  brevity,  we  adopt  a  loose  notation  of 
referencing  physical  components  by  subscripts. 

The  components  of  vector  (4)  define,  respectively,  conservation  of  mass,  conservation  of  momentum 
in  the  three  dimensions,  and  conservation  of  energy.  The  governing  system  is  closed  by  the  equation 
of  state.  For  this  work,  we  assume  that  the  fluid  is  air  and  behaves  as  a  perfect  gas,  whereby 

yMfp  —  pT,  (9) 

^  +  w^),  (10) 

7  —  1  2 

and  K  =  p/Pr. 

Four  dimensionless  parameters  emerge  from  nondimensionalization:  Mach  number,  Reynolds  num¬ 
ber,  Prandtl  number,  and  the  ratio  of  specific  heats,  defined,  respectively,  as  follows: 


Pr*Mr*^i*n 


Pr.W 

KT 


c* 


where  R*,  C*,  and  C*  are  the  ideal  gas  constant  and  the  specific  heats  at  constant  pressure  and 
constant  volume,  respectively.  Throughout  this  work,  dimensional  quantities  are  denoted  by  asterisk, 
reference  values  are  denoted  by  a  subscript  “r”,  Pr  =  0.7,  and  y  =  1.4.  In  the  nondimensionalization 
from  which  (4)  arises,  the  reference  values  for  density,  velocities,  and  temperature,  p*,  u*,  and  T* , 
respectively,  are  arbitrary.  Pressure  is  scaled  by  p*u*^.  Viscosity  is  normalized  by  the  viscosity  at  the 
reference  temperature  and  is  assumed  to  vary  according  to  Sutherland’s  law.  In  dimensionless  form 


+  C) 
T+C  ’ 


110.3  K 


Throughout  this  paper,  a  subscript  “e”  denotes  a  value  at  the  boundary-layer  edge,  and  lengths  are 
scaled  by  the  boundary-layer  displacement  thickness  at  the  inflow  boundary  5;*,  where,  in  general 
(White,  1974), 

I,  IR*)  Jo  V**/  \  pK)  Jl.  axisymmetnc. 

For  comparison  with  results  obtained  from  spatial  LST  and  the  PSE  method,  we  define  the  bound¬ 
ary-layer  length  scale  as  _ 


and  we  denote  the  Reynolds  number  based  on  edge  conditions  and  L*  as  Rbi^.  For  comparisons  with 
experiments,  it  is  also  useful  to  define  the  Reynolds  number  based  on  x*,  whereby  R^y.  =  {R^i)  • 
With  proper  interpretation,  (4)  is  valid  for  either  two-dimensional  or  axisymmetnc  bodies.  The 
two-dimensional  case  is  recovered  as  r  1  in  (3)  and  1/r  ^  0  in  H  (5). 
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A  nonconservative,  but  computationally  useful,  alternative  form  of  the  energy  equation  is 

Q4  =  p,  £4  =  pu,  £4  =  pv,  G4  =  pw, 

H4  =  a)-(r-i)pv*u,  (15) 

V-u  =  D^u  +  D^v  + 

where  repeated  indices  imply  summation.  In  subsequent  discussions  we  refer  to  (15)  as  the  “pressure 
equation”  as  distinguished  from  the  “energy  equation.” 

Mathematically,  DNS  is  the  numerical  solution  of  an  Initial-Boundary- Value  Problem  (IBVP),  For 
spatial  DNS,  inflow  and  outflow  boundary  conditions  are  required.  The  imposition  of  initial  and 
boundary  conditions  is  discussed  more  fully  in  the  next  section. 


3.  Numerical  Method 

Navier-Stokes  codes  fall  roughly  into  two  classes  depending  upon  the  application:  aerodynamic 
codes,  in  which  body  geometry  is  usually  complicated,  but  from  which  one  typically  wants  to  extract 
only  mean  quantities  such  as  surface  pressure,  lift,  and  drag;  and  DNS  codes,  in  which  body  geometry 
is  simple,  but  from  which  one  wants  to  compute  fine  details  of  the  flow  field,  ideally,  to  the 
smallest  scales.  Usually,  aerodynamic  codes  are  of  relatively  low-order  accuracy,  and,  to  capture 
shocks,  they  impose  significant  artificial  dissipation  through  some  form  of  upwinding.  Often  the 
boundary  layer  is  severely  underresolved,  a  deficiency  partially  atoned  for  by  the  use  of  transition  or 
turbulence  models.  In  contrast,  in  DNS,  for  which  the  boundary  layer  is  the  primary  focus,  dissipation 
and  dispersion  errors  must  be  minimized  if  the  growth  rates  and  phase  speeds  of  instability  waves  are 
to  be  computed  accurately. 

In  this  section,  we  present  a  fully  explicit  method  designed  specifically  for  the  DNS  of  instability 
and  transition  in  high-speed  boundary-layer  flows.  The  flow  region  of  interest  is  downstream  of  the 
bow  shock,  and  it  is  an  implicit  assumption  that  there  are  no  shocks  in  the  domain.  To  minimize 
dissipation  and  dispersion  errors,  we  rely  on  a  combination  of  spectral  methods  and  high-order 
central  compact-difference  approximations  for  spatial  derivatives.  In  the  development  of  the  algo¬ 
rithm,  accuracy,  efficiency,  and  simplicity,  in  that  order,  have  been  our  guiding  criteria. 

Spatial  Discretizations 

We  consider  the  physical  domain  and  grid  points  defined  by 

Xi„<Xi<Xout,  f  =  0, 

0<6»,<^,  jO,...,Ne,  (16) 

0<Zfc<z^3„  k  =  0, 

Let  Uijk  be  the  discrete  approximation  of  u{Xi,  6p  Zj^),  etc.,  and,  for  the  moment,  assume  that  grid  points 
are  equally  spaced  in  terms  of  their  respective  coordinates.  For  DNS  of  axisymmetric  (two-dimen¬ 
sional)  boundary  layers,  it  is  conventional  and  physically  reasonable  to  assume  that  the  flow  is 
periodic  in  the  azimuthal  (spanwise)  direction.  This  periodicity  permits  the  exploitation  of  spectral 
collocation  methods  (e.g.,  Canuto  et  al,  1988)  based  on  finite  series  expansions  of  the  flow  quantities 
in  terms  of  Fourier  basis  functions,  as,  for  example, 

Nell 

w(x,  0=  X  Mj.(x,z,  O-exp{/;n0},  (17) 

j-^-NeH 

where  i  =  Collectively,  the  Fourier  coefficients  Hj  provide  the  Fourier  spectrum,  an  intrinsic 

measure  of  the  adequacy  of  resolution.  The  integer  n,  which  appears  in  (16)  and  (17),  is  a  parameter  of 
the  flow  that  defines  the  period  in  the  azimuthal  direction  and  is  related  to  the  azimuthal  wave 
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number  jS  through  the  relation 

n  =  Pr  =  poRo,  Ro  =  R{xJ-  (18) 

(For  bodies  on  which  the  radius  changes,  such  as  a  cone,  (18)  implies  that  the  wave  number  of  oblique 
(helical)  waves  evolves  with  x.)  We  implement  the  collocation  method  in  the  conventional  pseudo- 
spectral  manner.  That  is,  derivatives  with  respect  to  6  of  the  truncated  series  in  (17)  are  evaluated 
exactly  in  Fourier  space,  whereas  nonlinear  terms  of  the  governing  equations  are  evaluated  in  physical 
space.  Vectorized  fast  Fourier  transforms  (FFTs)  are  used  to  shuttle  efficiently  between  the  transform 
and  physical  spaces.  If  desired,  aliasing  errors  are  controlled  by  spectral  truncation.  For  computa¬ 
tional  efficiency,  an  option  exists  in  the  code  to  enforce  symmetry  about  the  plane  0  =  0,  in  which 
case  odd  (e.g.,  v)  and  even  (e.g.,  u)  functions  are  expanded  in  sine  and  cosine  series,  respectively,  rather 
than  in  complex  exponential  series.  As  expected,  there  is  roughly  a  factor  of  2  reduction  in  computa¬ 
tional  effort  with  symmetry  enforced. 

The  proper  spectral  expansion  for  a  two-dimensional,  rather  than  an  axisymmetric,  body  is 
recovered  from  (17)  by  defining  the  azimuthal  arc  length  y  =  rd  and  by  using  (18). 

For  spatial  DNS,  both  the  streamwise  and  wall-normal  directions  are  aperiodic,  which  precludes 
the  use  of  Fourier  spectral  methods.  In  the  aperiodic  directions  our  DNS  code  allows  a  variety  of 
differentiation  options.  Among  these  are  a  Chebyshev  spectral-collocation  scheme  and  several  mem¬ 
bers  from  a  class  of  fourth-order  and  sixth-order  central  compact-difference  schemes  (Lele,  1990; 
Carpenter  et  al,  (1991)).  For  explicit  time  advancement,  the  Chebyshev  spectral  method  is  subject  to 
an  extremely  severe  restriction  on  the  time  step  of  the  form  At  ^  1/Nf,  which  renders  the  method 
impractical  for  long-time  integrations.  Nevertheless,  the  method  is  useful  for  diagnostic  purposes.  In 
particular,  it  permits  us  to  compute  the  Chebyshev  spectrum  as  an  indicator  of  the  adequacy 
of  resolution. 

In  the  context  of  fully  explicit  time  advancement,  derivatives  are  computed  by  means  of  compact- 
difference  techniques  as  follows: 

Mu^  =  £u,  (19) 

where  M  is  a  tridiagonal  matrix,  £  is  a  banded  matrix,  and  the  vector  u^,  for  example,  contains  the 
discrete  approximation  of  dujdz.  Matrices  M  and  £  are  referred  to  as  the  implicit  and  explicit 
operators,  respectively.  For  the  fourth-order  Fade  method,  £  is  tridiagonal.  For  sixth-order  methods, 
it  is  usually  pentadiagonal.  Compact-difference  operators  are  thus  of  the  form  C  =  M  and  are 
global  (dense  matrices). 

An  area  of  active  research  concerns  stable  boundary  closures  for  compact-difference  methods.  In 
general,  stencils  are  modified  in  the  vicinity  of  boundaries,  which  results  in  some  loss  of  formal 
accuracy.  It  is  well  known  (Gustafsson,  1975)  that,  for  a  hyperbolic  system  of  equations,  a  scheme  of 
hth  order  closure  can  retain  global  formal  accuracy  of  order  h  1  at  best.  In  distinguishing  between 
various  schemes,  we  adopt  the  nomenclature  of  Carpenter  et  al.  (1991).  For  example,  3,4-6-4,3  refers 
to  a  scheme  that  is  of  sixth-order  accuracy  at  interior  node  points,  fourth-order  accuracy  at  nodes  1 
and  JVj  —  1,  and  third-order  accuracy  at  nodes  0  and  N^.  Recently,  Carpenter  et  al.  (1991)  have 
developed  a  5,5-6-5,5  scheme  which  retains  global  sixth-order  accuracy.  It  evaluates  derivatives  at 
the  boundary  points  and  at  immediately  adjacent  points  by  fully  explicit  eight-point  stencils.  Among 
the  schemes  which  we  have  tried  and  for  which  options  exist  in  the  DNS  code  are  the  3—4-3  (Fade), 
3^4_6-4,3  (Lele),  and  5,5-6-5,5  (Carpenter  et  al.)  schemes.  The  reader  is  referred  to  Lele  (1990)  and 
Carpenter  et  al.  (1991)  for  details,  including  the  coefficients  of  matrices  M  and  £.  These  options  exist 
for  differentiation  in  both  the  streamwise  and  wall-normal  directions.  To  date  we  have  been  unable  to 
maintain  numerical  stability  whenever  fifth-order  boundary  closure  is  implemented  in  both  the  x  and 
the  z  directions.  Our  most  accurate  results  exploit  the  method  of  Carpenter  et  al.  in  the  streamwise 
direction  and  the  3,4-6-4,3  method  of  Lele  in  the  wall-normal  direction. 

All  differentiation  operators  assume  that  grid  points  are  equally  spaced  in  computational  space. 
Analytic  functions,  described  in  detail  later,  are  used  to  map  from  computational  space  onto  physical 
space,  where  the  grid  may  be  highly  stretched.  Metrics  of  the  transformation  are  incorporated  directly 
into  the  differentiation  operators. 
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Time  Advancement 

To  be  practical,  an  explicit  DNS  algorithm  must  provide  long-time  temporal  accuracy  while  requiring 
modest  temporary  storage.  Unless  excruciatingly  well  resolved  in  time,  we  believe  that  second-order 
schemes  lack  sufficient  accuracy  for  DNS.  Attempts  to  gain  higher-order  accuracy  through  multistep 
methods  violate  the  latter  practical  constraint.  An  elegant  compromise  is  afforded  by  the  use  of  one  of 
the  family  of  third-order  low-storage  Runge-Kutta  methods  proposed  by  Williamson  (1980).  To 
implement  the  Runge-Kutta  scheme,  we  cast  the  governing  equations  in  the  following  form: 

^  =  v 

dt 

(20) 

V  s  -D^E  -  D^F  -D^G-  H. 

Time  advancement  of  the  discrete  version  of  (20)  is  accomplished  in  physical  space.  Such  schemes  are 
now  in  wide  use,  and  so  we  omit  details. 

An  important  peripheral  issue  for  DNS  concerns  the  determination  of  a  limit  on  the  time  step  At, 
for  which  the  time-advancement  scheme  remains  stable,  but  which  is  not  unnecessarily  restrictive.  For 
parameter  values  typical  of  high-speed  boundary-layer  flows,  advection  and  viscous  constraints  on  the 
time  step  are  of  the  same  order  of  magnitude.  Moreover,  if  the  flow  is  transitional,  determination  of 
an  appropriate  time  step  is  complicated  by  localized  large-scale  velocity  and  temperature  fluctuations. 
Consequently,  we  have  chosen  to  estimate  the  stability  limit  dynamically  so  as  to  maintain  the  time 
step  continually  near  its  maximum  allowable  size.  Specifically,  we  require 

max(limi,  lim2)  *  At  <  sf, 

limi  =  (21) 

lim2  =  (lim^  + 

where  0  <  ^  <  1  is  the  stability  limit  safety  factor  (typically  sf  =  0.95),  and  where  lim^,  lim^,  and 
lim^  are  limiting  values  obtained  by  independent  consideration  of  the  stability  of  model  linear 
advection,  viscous  diffusion,  and  thermal  diffusion  equations,  respectively.  For  example,  to  evaluate 
the  advection  limit,  we  consider  the  linear  transport  equation 

<lt  +  (|w|  +  c)q^  +  (l^i  +  +  (|w|  +  c)q,  =  0.  (22) 

The  coefficients  u,  v,  and  w  are  obtained  from  the  instantaneously  “frozen”  velocity  fields  of  the  true 
problem.  Similarly,  the  speed  of  sound  c  —  y/ypip  is  computed  from  the  “frozen”  density  and  pressure 
of  the  true  problem.  The  integration  of  (22)  is  stable  provided 


At -lima  <  1, 


ijk  L 


'ijk  i\v\  +  c)ijk  (|w|  +  c)ijt 
i  cflyAyj  cfl.Az^  _ 


where  Ax^  =  —  x^.j,  etc.  The  Courant-Friedrichs-Levy  (CFL)  limits  cfl^,  cfly,  and  cfl^  depend  upon 

both  the  method  of  temporal  advancement  and  the  eigenvalue  spectrum  of  the  discrete  spatial 
operator  exploited  in  the  x,  y,  and  z  dimensions,  respectively.  Table  1  presents  the  stability  limits  of 
a  variety  of  spatial  operators  for  the  one-dimensional  linear  advection  (wave)  equation  integrated  in 
time  by  a  third-order  Runge-Kutta  (RK3)  scheme.  Some  of  these  limits  have  been  derived  by  rigorous 


Table  1.  Courant-Friedrichs-Levy  (CFL)  limit  |a(At/Ax)l  <  cfl  for 
the  scalar  linear  advection  equation  qt  +  aq^  =  0. 


Spatial  operator 

Stability  limit  cfl 

Method  of  derivation 

Compact  fourth  order 

1.0 

Empirical 

Compact  sixth  order 

^3/2 

Analytical 

Fourier  spectral 

Analytical 

Chebyshev  spectral 

Empirical 

Spatial  Direct  Numerical  Simulation  of  High-Speed  Boundary-Layer  Flows,  I 


59 


Table  2.  Viscous  stability  limit  |v(At/Ax^)|  <  vl  for  the  scalar  diffu¬ 
sion  equation 


Spatial  operator 

Stability  limit  vl 

Method  of  derivation 

Compact  fourth  order 

2.51/3 

Empirical 

Compact  sixth  order 

2.51/4 

Empirical 

Fourier  spectral 

2.51/7t^ 

Analytical 

Chebyshev  spectral 

2(2.51/n^) 

Empirical 

eigenvalue  analyses  and  are  so  noted.  Others  have  been  established  empirically  by  numerical  experi¬ 
mentation  with  the  time  step,  in  which  cases  the  scheme  “blows  up”  whenever  the  time  step  exceeds 
the  value  given.  For  the  Fourier  spectral  and  finite-difference  schemes,  the  stability  limits  are  sharp. 
For  the  Chebyshev  spectral  scheme,  the  limit  was  harder  to  evaluate  empirically  because  of  the  less 
explosive  nature  of  the  numerical  instability;  hence,  we  urge  some  caution  in  the  use  of  the  value 
presented. 

Similarly,  lim^  and  lim^  are  defined  by  considering  model  diffusion  equations  of  the  form 


Qt  = 


(24) 


with  the  help  of  Table  2.  The  domain  of  absolute  stability  for  third-order  Runge-Kutta  methods  can 
be  found  in  many  references  including  Canuto  et  al.  (1988).  The  factors  of  2.51  and  y/3  which  arise  in 
Tables  2  and  1,  respectively,  pertain  to  the  intercepts  of  the  RK3  stability  boundary  with  the  real  and 
imaginary  axes,  respectively.  The  stability  of  a  diffusion  problem  is  limited  by  the  value  —2,5\ 
of  the  real-axis  intercept,  whereas  the  stability  of  an  advection  problem  is  limited  by  the  values  ±^3 
of  the  imaginary-axis  intercepts.  We  conjecture  that  solutions  of  problems  of  mixed  advection- 
diffusion  type  remain  stable  provided,  say,  the  point  (  —  2.51  *lim^,  -lima)  lies  within  the  RK3 
domain  of  absolute  stability.  Equation  (21)  is  a  convenient  approximation  to  this  criterion  based  on 
the  fact  that  a  half-ellipse  through  the  real  and  imaginary  intercepts  lies  entirely  within  this  region. 

The  numerical  stability  of  semidiscrete  and  fully  discrete  compact-difference  schemes  is  subtle  and, 
as  mentioned  previously,  is  currently  an  area  of  active  research  (Carpenter  et  a/.,  1991).  In  general, 
what  is  true  for  scalar  equations  is  not  necessarily  true  for  systems.  Our  method  of  estimating  stability 
limits  is  somewhat  heuristic,  but  it  has  proved  a  useful  guideline  in  practical  applications  for  a  wide 
variety  of  numerical  methods  and  over  a  wide  range  of  parameter  values.  For  the  range  of  parameter 
values  of  the  test  cases  in  Section  5,  approximately  1000  time  steps  per  disturbance  period  suffice  to 
maintain  numerical  stability. 


Initial  Condition 

We  obtain  the  initial  (“basic”)  flow  state  Q^,  which  is  assumed  to  be  two-dimensional  or  ax- 
isymmetric,  from  the  spectrally  accurate  boundary-layer  (BL)  code  described  in  Pruett  and  Streett 
(1991)  and  Pruett  (1993).  Velocity  and  temperature  profiles  provided  by  this  code  are  smooth  nearly 
to  machine  precision.  Either  isothermal  or  adiabatic  wall  cases  can  be  computed  by  the  BL  code.  The 
present  results  assume  that  the  wall  is  adiabatic.  Because  the  boundary-layer  equations  neglect  terms 
of  order  and  higher,  a  small  residual  remains  when  the  boundary-layer  solution  is  injected  into 

the  steady  Navier— Stokes  equations.  Without  special  treatment,  this  numerical  error  introduces 
transients  into  the  domain  that  may  contaminate  the  Fourier  analysis  of  the  flow  field.  There  are  two 
philosophies  for  treating  the  initial  residual.  The  governing  equations  can  be  integrated  in  time 
(relaxed)  from  the  initial  guess  until  the  residual  becomes  arbitrarily  small  (in  practice,  on  the  order  of 
machine  zero).  Unfortunately,  for  an  explicit  code,  the  computation  of  a  clean  steady  state  may 
require  a  greater  computational  effort  than  the  computation  of  the  perturbed  flow  of  interest,  since 
numerical  errors  may  mimic  slowly  traveling  physical  waves.  Alternately,  small  source  terms  can 
be  subtracted  from  the  governing  equations  to  cancel  the  initial  residual  exactly  (Erlebacher  and 
Hussaini,  1990).  To  conserve  computational  resources,  we  adopt  the  latter  method.  For  the  parameter 
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subtracted  from  Navier-Stokes  equations  to  cancel  steady- 
Z  state  residual  of  base  flow. 


values  of  Cases  2  and  3  of  Section  5,  Figure  3  compares  the  magnitude  of  the  components  of  the 
source  term  f  that  are  subtracted  from  the  continuity,  u  and  w  momentum,  and  pressure  equations 
at  the  inflow  station.  The  source  term  is  small,  but  most  significant,  in  the  wall-normal  momentum 
equation,  where  the  slight  transverse  pressure  gradient  is  neglected  by  the  boundary-layer  equations. 
An  advantage  of  this  approach  is  that  it  facilitates  direct  comparisons  with  LST  and  PSE  analyses, 
both  of  which  assume  inherently  that  the  base  state  is  steady  and  simultaneously  satisfies  the 
governing  equations. 

For  the  presentation  of  results,  it  is  often  useful  to  display  only  the  disturbance  fields,  which  we 
denote  by  primes.  That  is,  Q'  =  Q  —  and  similarly  for  individual  components. 

Boundary  Conditions 

As  any  practitioner  of  spatial  DNS  will  acknowledge,  the  specification  of  boundary  conditions  for 
inflow- outflow  problems  is  a  delicate  matter.  The  reader  is  referred  to  Nordstrom  (1989)  and  Poinsot 
and  Lele  (1992)  for  detailed  discussions  of  the  issues  involved.  For  Navier-Stokes  calculations,  all 
flow  quantities  can  be  specified  at  the  inflow  boundary.  The  prescription  of  the  inflow  condition 
depends  upon  the  method  by  which  disturbances  are  introduced  into  the  flow.  To  date,  spatial  DNS 
calculations  have  introduced  forcing  either  through  a  time-periodic  inflow  condition  (e.g.,  originated, 
we  believe,  by  DeSanto  and  Keller,  (1962))  or  through  a  time-periodic  wall  boundary  condition  (e.g.. 
Krai,  1988;  Thumm,  1991).  For  stability  calculations,  the  difficulty  with  the  latter  method  is  that  it  is 
never  known  precisely  what  disturbances  are  generated.  Indeed,  in  high-speed  boundary-layer  flows,  it 
appears  that  multiple  modes  at  the  same  frequency  can  he  generated  by  periodic  suction  and  blowing 
at  the  wall  (Eissler  and  Bestek,  1993),  Therefore,  we  favor  the  former  approach.  At  the  inflow 
boundary  x  =  Xjn,  the  flow  is  specified  as  the  superposition  of  the  steady  two-dimensional  or 
axisymmetric  base  flow  and  a  temporally  periodic  fluctuation  of  frequency  o,  amplitude  e,  streamwise 
wave  number  a,  and  spanwise  (azimuthal)  wave  number  P  as  follows: 

Q(Xi„,  y,  z,  t)  =  z)  +  a  exp[i(aXi„  +  Py)  -  (of]'P(z)  +  c.c.  (25) 

To  minimize  temporal  transients,  £  is  ramped  smoothly  to  its  maximum  value  over  an  interval  of  time 
(typically  one  disturbance  period).  The  disturbance  structure,  which  is  contained  in  the  complex  vector 
T,  is  obtained  either  from  spatial  LST  (Ng  and  Zang,  1993)  or  from  PSE  theory  as  in  Chang  et  al 
(1991).  In  either  case  the  disturbance  is  normalized  so  that  the  maximum  amplitude  of  the  tempera¬ 
ture  fluctuation  is  unity.  For  high-speed  flows,  the  disturbances  derived  from  the  parallel  (LST)  and 
nonparallel  (PSE)  theories  are  significantly  different,  as  shown  in  Figure  4,  and  the  disparity  tends  to 
increase  with  x  and  with  Mach  number.  For  developing  (nonparallel)  high-speed  boundary-layer 
flows,  the  use  of  LST-derived  forcing  results  in  a  significant  inconsistency  at  the  inflow  boundary, 
which  manifests  itself  in  undesirable  streamwise  transients  downstream  of  the  inflow  boundary. 
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Figure  4.  Amplitude  of  temperature  component  of  linear  dis¬ 
turbance  at  ^  700  from  parallel  (LST)  and  nonparallel 
(PSE)  theories. 


Re^ 

Figure  5.  Transients  downstream  of  inflow  boundary  for  Case 
4  of  Section  5. 


Moreover,  if  the  amplitude  of  the  disturbance  is  large  enough,  a  weak  Mach  wave  can  be  seen 
emanating  from  the  inflow  boundary  near  the  critical  layer.  It  is  therefore  important  that  nonparallel 
effects  be  considered.  We  have  found  PSE-derived  forcing  functions  virtually  to  eliminate  inflow 
inconsistencies,  as  shown  in  Figure  5,  which  presents  the  evolution  of  the  temperature  maximum  in 
the  region  near  the  inflow  boundary  for  Case  4  in  Section  5.  In  contrast,  the  wandering  of  the  tem¬ 
perature  maximum  away  from  the  nearly  constant-slope  DNS/PSE  curve  for  the  DNS/LST  case  is  symp¬ 
tomatic  of  the  inconsistency  between  the  parallel  LST  theory  and  the  nonparallel  DNS  calculation. 

For  reference,  the  amplitudes  of  each  of  the  five  components  'P  of  an  oblique  second-mode 
disturbance  in  a  Mach  6.8  boundary-layer  flow  on  a  cone  are  shown  in  Figure  1,  having  been 
obtained  from  the  nonparallel  PSE  method  for  parameter  values  which  correspond  to  Case  4  in 
Section  5.  Note  the  domination  of  the  temperature  and  density  components  of  the  disturbance. 

In  the  vicinity  of  the  outflow  boundary,  we  have  found  a  buffer  domain,  proposed  by  Streett  and 
Macaraeg  (1989/90),  to  be  effective  in  passing  large-amplitude  ff actuations  with  minimal  reflection  and 
upstream  influence.  Within  the  buffer  region,  the  governing  system  of  equations  is  gradually  modified 
as  follows: 

1.  Streamwise  viscous  terms  are  smoothly  attenuated  to  zero  to  parabolize  the  governing  system. 

2.  The  base  streamwise  velocity  profile  is  smoothly  brought  to  that  of  a  uniform  flow  at  the 
velocity  of  the  free  stream  to  ensure  that  all  characteristics  lead  out  of  the  domain. 

3.  The  source  term  f  discussed  previously  is  modified  significantly  in  the  buffer  region  to  balance 
the  changes  to  the  basic  flow. 

With  these  changes,  flow  quantities  along  the  outflow  boundary  can  be  extrapolated  from  the  interior. 

At  the  wall,  the  usual  no-slip  conditions  are  imposed  on  the  velocities  (except  within  the  buffer 
region).  Temperature  is  assumed  to  be  fixed  at  its  adiabatic  value.  (The  rationale  for  the  hybridization 
of  the  adiabatic  and  Dirichlet  wall  condition  on  temperature  is  addressed  in  an  article  by  Pruett  and 
Zang  (1992)  in  a  previous  issue  of  this  journal.)  No  condition  is  imposed  on  density;  its  value 
is  determined  directly  from  the  governing  equations  at  the  wall.  Pressure  at  the  wall  is  derived  from 
density  and  temperature  via  the  equation  of  state. 

At  the  far-field  boundary  (z  =  we  adapt  the  nonreflecting  boundary  conditions  proposed  by 
Thompson  (1987),  which  are  based  in  inviscid  characteristic  theory.  Here,  these  are  applied  only  to  the 
disturbance  field  U'  =  U  -  in  the  following  way.  Adding  and  subtracting  formally  identical 
quantities  from  the  right-hand  side  of  (20),  we  obtain 

dQ  dG'  5G 

dt  dz  dz 


(26) 
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where 

dG'  _  dG  dU 
dz  d\J  dz  * 

Recognizing  that 

5Q_5Q5U 


(27) 


and  exploiting  the  notation  A  =  5Q/5U  and  B  =  dG/dV,  we  temporarily  recast  (26)  in  terms  of  U  as 
follows: 


dV 

dt 


dz 


(28) 


where  S  =  A~^B,  and  S  has  yet  to  be  defined.  Matrix  S  is  then  diagonalized  by  the  similarity 
transformation  S  =  PAP"\  where  A  is  the  diagonal  matrix  with  eigenvalues  {w  —  c,  w,  w,  w, 
w  +  c}  in  the  order  shown.  Positive  and  negative  eigenvalues  correspond  to  outbound  and  inbound 
characteristics,  respectively.  Now  let  S  =  PAP"\  where  A  is  the  diagonal  matrix  with  eigenvalues 
=  max(A^,  0).  Writing  (28)  once  again  in  terms  of  Q,  we  obtain  our  working  form 


^  =  V-^P[A-A]P”^ 

Note  that  if  U  =  U®,  or  if  all  eigenvalues  of  A  are  nonnegative,  then  there  is  no  change  to  the 
right-hand  side  vector  V.  If,  however,  there  are  disturbances  traveling  inbound  at  the  far-field 
boundary,  their  time  derivatives  are  set  to  zero  on  the  boundary. 


Mappings 

As  shown  in  Figure  1,  the  eigenfunction  of  a  second-mode  disturbance  in  a  high-speed  boundary-layer 
flow  has  a  double  structure  in  which  the  temperature  component  is  sharply  peaked  near  the  wall  and 
the  critical  layer.  Consequently,  accuracy  and  resolution  considerations  necessitate  that  grid  points  be 
clustered  in  both  regions  of  sharp  gradients.  For  this  purpose,  we  use  a  highly  tuned  mapping 
from  the  computational  space  to  the  physical  space.  The  mapping  is  adapted  from  Erlebacher  and 
Hussaini  (1990)  and  combines  gradual  exponential  stretching  away  from  the  wall  with  clustering  in 
the  vicinity  of  the  critical  layer.  In  computational  space,  for  fc  =  0,  1, ...,  N^, 


compact  difference, 
Chebyshev  spectral. 


(30) 


For  the  compact-difference  schemes,  note  that  grid  points  are  equally  spaced  in  computational  space. 
There  are  five  parameters  for  the  mapping:  Zy2^  Zq,  Az,  and  t.  The  following  transcendental 

equation  maps  the  interval  [—1,  +1]  in  computational  space  onto  itself  and  performs  a  clustering  of 
grid  points  of  strength  t  and  width  Az  about  the  point  Zq: 

C  +  Ttanh  -1<^<+1.  (31) 

I  J 


The  subsequent  exponential  transformation  below  maps  the  computational  interval  [—1,  +1]  onto 
the  physical  interval  [0,  z^^^y. 
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Figure  6.  Wall-normal  (a)  mapping  and  (b)  associated  metric. 


The  four  remaining  unknown  quantities  in  (31)  are  determined  by  requiring 

z(Co)  =  ^0.  C(-l)=-l>  C(+l)=+l-  (33) 

If  T  =  0,  there  is  no  clustering  about  Zq,  in  which  case  exactly  half  of  the  grid  points  lie  between  the 
wall  and  Zy2-  We  find  (32)  to  be  superior  to  the  bilinear  fractional  transformation  used  originally  by 
Erlebacher  and  Hussaini  (1990),  which  stretches  too  fast  in  the  far  field.  Figure  6  shows  the  grid-point 
distribution  and  metric  for  the  parameter  values  of  Mach  4.5  in  Cases  2  and  3  in  Section  5.  In 
practice,  we  tune  the  number  of  grid  points  and  mapping  parameters  by  use  of  the  temporal  DNS 
code  of  Pruett  and  Zang  (1992).  More  specifically,  we  adjust  the  resolution  and  mapping  parameters 
until  we  are  able  to  recover  local  (global)  disturbance  growth  rates  from  the  DNS  that  are  in  four 
(six)  place  agreement  with  eigenvalues  obtained  from  temporal  LST. 

Similarly,  a  streamwise  mapping  can  be  used  to  concentrate  points  downstream  where  nonlinear 
interactions  lead  to  a  broadening  of  the  spectrum  in  wave  space.  For  this  paper,  however,  all  results 
were  obtained  with  grid  points  equally  spaced  in  x. 

In  spatial  DNS  the  boundary  layer  thickens  as  x  increases.  Eventually,  the  mapping  described 
above  will  become  detuned,  and  the  far-field  boundary  at  will  pinch  the  flow  unless  the  physical 
domain  grows  in  wall-normal  extent  approximately  as  fast  as  the  boundary  layer.  Motivated  by 
boundary-layer  theory,  we  define  z  =  flf{x)  where  fix)  is  a  smooth  function  which  grows  like  the 
boundary  layer.  For  laminar  flow,  an  appropriate  choice  is  fix)  =  ^x/xq.  An  option  is  provided  in 
the  DNS  code  to  fix  the  far-field  extent  of  either  z  or  r\.  In  the  latter  case  the  physical  grid  is 
nonorthogonal  and  streamwise  and  wall-normal  derivatives  are  modified  via  the  chain  rule  as  follows: 

dx'^dx  fix)  dtj’ 
d  1  d 
dz  fix)dri' 

With  the  exception  of  the  first  case  in  Section  5,  the  results  presented  below  were  obtained  with  the 
growing-domain  option. 


Filtering 

From  the  work  of  Trefethen  (1982),  Vichnevetsky  (1986),  and  Poinsot  and  Lele  (1992),  it  is  clear 
that  all  finite-difference  schemes  reflect  some  energy  at  outflow  boundaries.  Typically  the  reflection  co- 
efficient  (i.e.,  the  ratio  of  incident  to  reflected  energy)  behaves  like  (aAx)^  where  h  is  the  order  of 
the  scheme  and  a  is  the  disturbance  wave  number.  Moreover,  numerical  reflections  arise  even  if  the 
physical  boundary  conditions  imposed  are  perfectly  nonreflecting.  Regardless  of  the  order  of  the 
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Figure  7.  Two-point  oscillations  and  spurious  physical  waves  in  disturbance  pressure  p'  in  absence  of  filtering.  Two-dimen¬ 
sional  Mach  4.5  flat-plate  boundary-layer  flow  with  parameter  values  similar  to  Case  2  of  Section  5.  Streamwise  extent  limited 
to  nine  disturbance  wavelengths.  Wall-normal  extent  0  <  ^  <  7.5.  Upper  and  lower  photographs,  respectively,  show  instants  in 
time  before  and  after  exit  of  leading  wavefront  from  computational  domain. 


scheme,  energy  at  the  “sawtooth”  wavelength  (aAx  =  tt)  is  usually  totally  reflected.  As  shown  by 
Trefethen  (1982)  and  Vichnevetsky  (1986),  reflected  energy  travels  upstream  as  a  wave  packet  at  the 
group  velocity  which  depends  on  the  spatial  discretization  scheme.  As  mentioned  by  Poinsot  and 
Lele  (1992),  for  the  one-dimensional  Euler  equation,  the  group  velocities  of  the  fourth-order  and 
sixth-order  central  compact-difTerence  schemes  are  3(|u|  +  c)  and  14(|u| -f- c)/3,  respectively.  Once 
the  wave  packet  encounters  the  Dirichlet  inflow  condition,  which  is  numerically  reflecting  to  all 
wavelengths,  it  reforms  and  travels  downstream  as  a  spurious  physical  disturbance.  For  the  sixth- 
order  compact-difference  scheme,  this  entire  phenomenon  masquerades  as  an  apparent  coupling 
between  the  outflow  and  inflow  boundaries,  whereby  spurious  physical  oscillations  emanate  from  the 
inflow  boundary  shortly  following  the  exit  of  the  leading  wavefront  from  the  outflow  boundary,  as 
shown  by  the  flow-visualization  sequence  in  Figure  7.  The  photograph  shows  the  instantaneous 
disturbance  pressure  p'  at  two  different  times  in  a  simulation  of  a  two-dimensional  Mach  4.5  flat-plate 
boundary-layer  flow.  The  parameter  values  are  similar  to  those  of  Case  2  of  Section  5.  However,  the 
corrective  filtering  described  below  has  not  been  implemented,  and,  for  clarity,  the  domain  is 
foreshortened  to  a  length  of  nine  disturbance  wavelengths  in  streamwise  extent.  The  upper  photo¬ 
graph  shows  p'  prior  to  the  arrival  of  the  leading  wavefront  at  the  outflow  plane.  Two-point 
oscillations,  predominantly  in  the  wall-normal  direction,  are  clearly  visible  just  outside  the  boundary 
layer.  The  lower  photograph,  with  precisely  the  same  color  scale,  shows  p'  several  periods  after  the 
exit  of  the  leading  wavefront  from  the  domain.  Note  the  appearance  of  spurious  physical  waves  along 
the  upper  boundary,  the  ultimate  result  of  a  numerical  double  reflection  off  the  outflow  and  inflow 
boundaries. 

One  of  several  possible  solutions,  including  upwind  biasing,  is  through  filtering,  which  we  prefer  for 
its  simplicity  and  efficiency.  Outside  the  boundary  layer,  the  physical  viscosity  has  a  vanishingly  small 
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Figure  8.  Transfer  function  tr(aAx)  for  sixth-order  compact- 
difference  filter. 

effect;  without  some  added  numerical  dissipation,  the  central-difference  approximations  result  in 
odd-even  decoupling  that  introduces  energy  at  the  troublesome  “sawtooth”  wavelength.  The  trick  is 
to  keep  energy  at  this  wavelength  from  crossing  the  outflow  boundary.  Here,  minimal  damping  is 
imposed  by  periodically  applying  a  low-pass  sixth-order  compact-difference  filter  (Lele,  1990)  to  the 
solution.  Typically  we  apply  the  filter  every  few  time  steps,  in  which  case  the  additional  computational 
effort  is  insignificant  when  amortized  over,  say,  nine  Runge-Kutta  stages.  Moreover,  we  find  it 
necessary  to  apply  the  filter  only  in  the  wall-normal  direction. 

Lele’s  sixth-order  filter  has  two  free  parameters.  For  the  values  we  have  chosen,  the  formal 
truncation  error  for  filtering  a  function  f{x)  is  and  the  transfer  function  is  shown  in 

Figure  8.  Two-point  oscillations  (shown  by  the  vertical  dotted  line)  are  completely  eliminated  by  the 
filter,  whereas  oscillations  of  twice  that  wavelength  (shown  by  the  vertical  dashed-dotted  line)  and 
longer  are  virtually  undamped.  In  practice,  the  filter  incorporates  third-order  boundary  closure  (Lele, 
1990),  which  reduces  the  global  accuracy  to  fourth  order,  as  does  the  boundary  closure  of  the 
3^4_6_4^3  compact-difference  scheme.  The  filter  parameters  and  frequencies  in  current  use  are  in  no 
sense  optimal,  and  it  may  well  be  possible  to  reduce  the  modest  computational  effort  still  further  by 
less  frequent  filtering. 

An  additional  subtlety  arises  when  implementing  the  filter  in  the  context  of  a  total  variable 
formulation  of  the  governing  equations  (rather  than  a  disturbance  variable  formulation).  Over  many 
applications,  filtering  effects  gradual  evolution  of  the  base  state.  For  example,  for  a  simulation  similar 
to  that  of  Case  4  in  Section  5,  but  with  a  foreshortened  domain  (one-sixth  the  streamwise  extent  of 
Case  4)  and  no  imposed  disturbance  (^  =  0,0),  the  base  density  changes  approximately  0.01%  in  1.4 
flow-through  times  when  filtering  is  applied  every  third  time  step.  This  presents  two  problems.  First, 
because  of  the  extreme  sensitivity  of  the  stability  of  the  flow  to  changes  in  the  base  state,  unintended 
evolution  of  the  base  state  may  affect  hydrodynamic  stability.  Second,  when  the  basic  state  is 
subtracted  from  the  total-flow  variables  for  purposes  of  flow  visualization,  the  disturbance  fields, 
whose  velocity  components  may  be,  say,  of  order  10“^,  are  contaminated.  The  difficulty  is  corrected 
by  incorporating  an  additional  steady  source  term  fi,  to  be  defined  shortly,  into  the  right-hand  side  of 
the  governing  equations.  For  simplicity,  in  the  context  of  a  first-order  (forward  Euler)  time-integration 
scheme,  the  filtering  algorithm  can  be  summarized  as  follows: 

For/-0,  1,  2,... 


I  +  At  •  [v(qO  -  f  -  fi]  (35) 

I  q^^^=Fq*^^ 

where  q  and  v  denote  the  discrete  representations  of  Q  and  V,  respectively,  a  tilde  denotes  a  filtered 
quantity,  F  is  the  discrete  filter  operator,  and  I  is  the  time-level  index.  In  practice,  Pq*  =  Hq\  where  P 
and  H  are  pentadiagonal  and  heptadiagonal  matrices,  respectively,  whereby  F  =  P~^H.  Recognizing 
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Figure  9.  Instantaneous  disturbance  pressure  p'  at  time  of  53  periods  of  oscillation  of  fundamental  for  Case  4  of  Section  5.  For 
clarity,  computational  domain  truncated  in  streamwise  extent  to  1436  <  <  1536.  (Buffer  domain  shown.)  Wall-normal 

extent  0  <  rj  <  7.5.  Results  shown  in  “peak”  plane  6^0.  Colors  depict  contour  levels  between  —  3  x  10“^  and  3  x  10“^, 


that  (by  definition)  v(q°)  —  f  =  0  identically,  the  reader  can  readily  show  from  (35)  that  if  /  =  0  and 
ft  =  [(F  —  /)/At]q°,  then  q^  =  q^.  That  is,  formally  the  base  state  is  preserved  over  the  first  and  all 
subsequent  time  steps  in  the  absence  of  forced  disturbances.  In  practice,  the  additional  source  term 
is  “turned  on”  only  immediately  prior  to  filtering  operations,  in  which  case  the  filtered  basic  state 
remains  constant  over  time  nearly  to  machine  precision. 

Figure  9  is  presented  to  illustrate  the  beneficial  effects  on  the  solution  of  the  filtering  algorithm 
described  above,  including  the  adjustment  to  prevent  unintentional  evolution  of  the  base  state.  It 
portrays  the  instantaneous  disturbance  pressure  p'  obtained  from  Case  4  of  Section  5.  For  clarity, 
however,  only  the  final  ten  disturbance  wavelengths  of  the  computational  domain  (including  the  buffer 
domain)  are  shown.  The  time  corresponds  to  53  periods  of  oscillation  of  the  disturbance.  We  note 
that  the  maximum  amplitude  of  p'  is  quite  small,  in  our  normalization,  more  than  two  orders  of 
magnitude  below  that  of  T\  whose  amplitude  at  the  inflow  boundary  is  e  =  0.001.  In  contrast  to 
Figure  7,  two-point  oscillations  have  been  completely  damped,  thereby  eliminating  detectable  bound¬ 
ary  reflections  and  spurious  physical  waves,  without  detriment  to  the  evolution  of  the  instability  wave. 

4.  Algorithm  Details 

In  algorithm  design  there  is  frequently  a  tradeoff  between  computational  effort  and  storage  require¬ 
ments,  as  is  the  case  here.  For  a  fully  explicit  scheme,  storage  requirements  are  modest,  and  the 
overriding  consideration  is  to  minimize  the  computational  work  per  time  step.  We  discuss  very  briefly 
several  considerations  which  contribute  to  the  efficiency  and  simplicity  of  the  algorithm. 

For  the  three-dimensional  compressible  Navier-Stokes  equations,  the  minimum  number  of  partial 
derivative  computations  required  by  an  explicit  algorithm  to  evaluate  V  is  27:  9  for  the  rate-of- 
deformation  tensor,  3  for  the  heat  flux  components,  and  the  remaining  15  for  the  components  of  the 
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flux  vectors  E,  F,  and  G.  Both  our  conservative  (energy  equation)  and  nonconservative  (pressure 
equation)  formulations  require  only  the  minimum  number  of  27  partial  derivative  evaluations.  Specifi¬ 
cally,  partial  derivative  evaluations  are  implemented  through  calls  to  a  subroutine  PARTIAL  whose 
arguments  specify: 

1.  The  direction  of  the  derivative. 

2.  Whether  or  not  to  include  certain  terms  which  arise  from  nonzero  Christoffel  symbols  as  in  the 
distinction  between  D^u  and  D^u  in  (3). 

3.  If  spanwise  (azimuthal)  symmetry  is  enforced,  whether  or  not  the  function  is  even  or  odd. 

Parameters  passed  to  PARTIAL  also  specify  the  difference  schemes  for  each  coordinate  direction.  By 
this  construction  we  have  a  wide  range  of  possible  options  without  complicating  the  core  subroutine 
that  evaluates  the  right-hand  side  vector  V.  The  present  algorithm  uses  (we  believe)  the  minimum 
storage  possible  given  the  minimum  work  constraint  previously  discussed.  The  conservative  (energy) 
formulation  of  the  governing  equations  requires  18  storage  arrays  of  approximate  size  A/^  x  x  N^, 
6  for  the  primitive  variables  U  and  the  viscosity,  5  for  storing  components  of  the  right-hand  side 
vector  V,  6  for  temporary  use  in  the  evaluation  of  V,  and  1  for  partial  derivative  evaluations.  A 
significant  advantage  of  the  nonconservative  (pressure)  formulation  is  that  only  16  storage  arrays  are 
required.  Since  the  energy  equation  and  pressure  equation  formulations  yield  essentially  identical 
results,  we  favor  the  pressure  formulation  for  computational  efficiency.  For  both  formulations,  one 
additional  array  is  required  if  the  growing-domain  option  is  invoked. 

Finally,  for  computational  efficiency,  innermost  loops  typically  range  over  all  {N^  +  1)  x  (Nq  +  1) 
grid  points  in  a  surface  of  constant  /c,  ensuring  nearly  optimal  speed-up  on  vector  processors.  For 
directionally  dependent  calculations  such  as  FFT  evaluations  and  the  solution  of  tridiagonal  systems, 
the  inner  computational  loop  ranges  over  the  number  of  independent  transforms  or  right-hand  sides. 
On  single  processors  of  the  Cray  Y-MP  and  the  Cray  C90,  the  present  algorithm  performs  at  greater 
than  150  megaflops  and  415  megaflops,  respectively.  In  recent  calculations  on  the  C90,  the  algorithm 
has  been  assessed  at  5.6  x  10“^  seconds  per  gridpoint  per  (full)  time  step. 


5.  Validation 

For  comparisons  with  DNS  results,  we  rely  on  compressible  LST  (Mack,  1984,  Ng  and  Zang,  1993) 
and  on  results  obtained  by  the  compressible  PSE  method  of  Chang  et  al.  (1991).  The  reader  is  referred 
to  these  authors  for  details.  Briefly,  in  classical  LST,  an  eigenvalue  problem  results  from  certain 
assumptions  about  the  basic  flow  and  the  wave-like  nature  of  the  disturbances.  Because  of  the 
assumptions,  the  results  of  LST  are  strictly  valid  only  for  low-amplitude  (linear)  disturbances  in 
parallel  base  flows.  Our  LST  results  were  obtained  from  the  spectrally  accurate  spatial  linear  stability 
code  of  Ng  and  Zang  (1993). 

The  PSE  approach  uses  a  traveling  wave  ansatz  similar  to  LST,  except  that  the  disturbance  shape 
function  (T'  in  (25))  is  allowed  to  vary  in  both  the  wall-normal  and  the  streamwise  directions.  Rapid 
oscillation  of  the  wave  is  incorporated  into  the  exponential  part  of  the  wave  ansatz,  whereby  the 
shape  function  evolves  in  x  on  a  scale  much  longer  than  a  wavelength.  The  governing  equations 
thereby  reduce  to  a  set  of  PDEs  for  the  shape  function  only.  The  PSE  method  becomes  approximate, 
rather  than  exact,  when  these  PDEs  are  parabolized  to  facilitate  a  marching  solution.  Provided  the 
instabilities  are  of  convective  nature,  as  they  are  for  most  high-speed  boundary-layer  flows,  the 
parabolization  approximation  is  quite  reasonable.  For  nonlinear  problems,  the  disturbances  are 
expressed  as  Fourier  series  in  the  frequency  domain.  The  equations  for  each  Fourier  mode  are 
independent  except  through  inhomogeneous  terms  which  arise  due  to  nonlinear  mode  interactions. 
Therefore,  the  PSE  method  can  treat  both  nonparallel  and  moderately  nonlinear  effects. 

The  comparison  of  DNS  results  with  those  obtained  by  LST  and  PSE  requires  the  Fourier 
decomposition  of  the  DNS  data  in  time  and  in  the  azimuthal  (spanwise)  dimension.  The  resulting 
modes  are  identified  with  ordered  pairs  m^),  where  the  integers  and  label  the  harmonics  in 
Fourier  space  with  respect  to  the  temporal  frequency  co  and  the  azimuthal  wave  number  respec¬ 
tively.  For  example,  the  fundamental  mode  is  labeled  (1,  0)  if  two-dimensional  or  (1, 1)  if  oblique. 
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R  *  in-2  Figure  10.  First-  and  second-mode  unstable  regions  for  Mach 

®L  ^  4.5  flat-plate  boundary  layer  (after  Mack,  1984). 


We  present  four  validation  cases.  The  first  three  consider  two-dimensional  second-model  distur¬ 
bances  in  Mach  4.5  flow  along  a  flat  plat.  Of  these  three,  the  first  two  examine  the  evolution  of 
low-amplitude  disturbances  in  parallel  and  nonparallel  boundary-layer  flows,  respectively.  In  the  third 
test  case,  forcing  is  of  large  amplitude  and  harmonics  generated  by  nonlinear  interactions  attain 
significant  amplitudes.  The  final  test  case  considers  the  evolution  of  a  pair  of  symmetric  oblique 
(helical)  second-mode  disturbances  in  a  Mach  6.8  boundary-layer  flow  along  a  cone. 

As  a  point  of  reference  for  the  flat-plate  cases,  we  include  Figure  10  (adapted  from  Mack,  1984), 
which  shows  distinct  regions  of  instability  for  first-  and  second-mode  disturbances  in  a  Mach  4.5 
planar  boundary-layer  flow.  However,  we  caution  against  attempts  at  exact  comparisons  for  these 
reasons:  Mack’s  diagram  is  derived  from  temporal,  rather  than  spatial,  LST;  LST  inherently  assumes 
locally  parallel  flow;  and  Mack  uses  a  slightly  different  formula  for  viscosity. 

Case  1:  Flat  Plate,  Mach  4.5  Parallel  Flow,  Two-Dimensional  Linear  Disturbance.  Our  primary  purpose 
here  is  to  evaluate  the  resolution  required  in  the  streamwise  direction  to  capture  accurately  the 
evolution  of  a  monochromatic  disturbance.  Figure  11  compares  the  computed  (DNS)  and  theoretical 
(LST)  maximum  amplitude  of  the  temperature  fluctuation  for  a  disturbance  of  dimensionless  fre¬ 
quency  F  =  cu*L*/(Me  =  2.29  X  10”\  Re^  ~  955,67,  and  amplitude  e  «  1  in  flat-plate  adiabatic- 
wall  flow  with  Mg  =  4.5  and  T*  =  61.11  K.  These  parameter  values  define  a  slightly  damped  second¬ 
mode  disturbance  that  corresponds  to  a  point  just  beyond  the  upper  branch  neutral  curve  in  Figure 
10.  For  the  DNS,  the  computational  domain  spans  eight  wavelengths  in  streamwise  extent  (based  on 
parallel,  linear  theory).  The  amplitude  maxima  versus  x  are  determined  by  Fourier  time-series  analysis 
over  the  interval  from  13  to  14  periods  of  oscillation  of  the  disturbance,  after  the  leading  wavefront 
has  exited  the  domain.  DNS  results  are  shown  for  the  5,5-6-5,5  scheme  at  streamwise  resolutions  of 
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Figure  11.  Case  1:  DNS  and  LST  results  for  linear  stability 
problem.  In  the  legend  “DNS  16/6,”  for  example,  denotes 
results  obtained  from  spatial  DNS  by  6th-order  differencing 
and  16  gridpoints  per  wavelength  in  the  streamwise  direction. 
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8,  12,  and  16  grid  points  per  disturbance,  wavelength  and  for  the  3-4-3  scheme  with  16  grid  points 
per  wavelength.  (In  the  legend  for  Figure  11,  the  first  and  second  integers  refer  to  the  number  of  grid 
points  per  wavelength  and  the  order  of  the  scheme,  respectively.)  In  all  cases,  wall-normal  differencing 
is  accomplished  by  the  3,4-6-4,3  compact-difference  scheme  with  =  144.  For  direct  comparison 
with  LSI,  the  base  flow  is  parallel  (i.e.,  u  =  u(z),  T  =  T{z),  and  w  =  0),  in  which  case  the  forcing  terms 
analogous  to  Figure  3  are  considerably  larger  than  in  the  developing  (nonparallel)  flow  case.  For  the 
fourth-order  scheme  and  for  the  sixth-order  scheme  with  fewer  than  12  points  per  wavelength, 
significant  oscillations  are  seen  in  the  maxima.  For  the  sixth-order  scheme,  the  eigenvalue  extracted 
from  the  DNS  is  correct  to  about  three  and  four  places  for  12-  and  16-point  resolutions,  respeetively, 
except  in  the  immediate  vicinity  of  the  inflow  boundary  where  there  is  a  very  slight  streamwise 
transient.  The  resolution  needed  for  accurate  results  is  about  twice  the  six  to  eight  points  per 
wavelength  we  would  have  naively  anticipated.  For  nonlinear  problems,  these  results  suggest  that  it  is 
desirable  to  resolve  the  shortest  wavelengths  (highest  harmonics)  with  at  least  12  points  per  wave¬ 
length. 

Figure  11  also  offers  reasonable  validation  of  the  buffer-domain  outflow  condition;  no  appreciable 
influence  is  evident  on  the  decay  rate  of  the  instability  wave  upstream  of  the  edge  of  the  buffer 
domain  (shown  by  the  vertical  line). 

Case  2:  Flat  Plate,  Mach  4.5  Nonparallel  Flow,  Two-Dimensional  Linear  Disturbance.  Case  2  corre¬ 
sponds  approximately  to  the  ray  that  extends  through  points  1  and  2  in  Figure  10  (Mack,  1984).  In 
this  case  an  evolving  (nonparallel)  Mach  4.5  flat-plate  boundary-layer  flow  with  the  same  flow 
parameters  as  before  is  forced  at  a  dimensionless  frequency  F  =  2.2  x  10”^.  We  note  that,  for  Cases 
1-3,  the  ratio  5*/L*  is  constant  along  the  plate  with  the  value  10.46.  For  the  DNS,  the  inflow 
and  outflow  boundaries  correspond  to  Re^  =  0.49  x  10®  {Rej^  =  700)  and  Re^  =  1.404  x  10®  {Rej^  = 
1185.3),  respectively.  Over  this  range  of  Reynolds  numbers,  the  selected  frequency  excites  a  second¬ 
mode  disturbance.  Relative  to  Figure  10,  the  computational  domain  is  36  wavelengths  long  (based  on 
linear,  nonparallel  theory  at  Xi„)  and  spans  roughly  from  point  1  (just  prior  to  the  lower  branch 
neutral  point)  to  point  2  (somewhat  beyond  the  upper  branch  neutral  point).  For  the  DNS,  =  144 
and  =  7.5.  Grid-point  clustering  and  stretching  is  such  that  75%  of  all  points  fall  in  the  region 
bounded  by  the  wall  and  z  =  1.4.  The  buffer  region  spans  the  last  5%  of  the  domain  in  x  (i.e.,  the  last 
1.8  wavelengths).  For  comparison  with  LST,  a  small  amplitude  of  e  =  0.001  is  used  for  the  DNS 
calculation  to  render  insignificant  the  nonlinear  terms  of  second  order  or  higher  in  e.  The  present 
results  were  obtained  for  =  576,  with  points  equally  spaced  in  x,  and  can  be  considered  very  well 
resolved  at  16  grid  points  per  wavelength. 

The  PSE  and  DNS  calculations  both  proceed  from  identical  disturbance  states.  However,  their 
respective  base  states  were  derived  independently,  that  of  the  PSE  calculation  being  obtained  by 
finite-difference  techniques.  The  amplitude  distributions  of  the  disturbance  components  imposed  at  the 
inflow  station  can  be  found  in  Figure  2  of  Pruett  and  Chang  (1993).  Because  the  PSE  calculation  is 
linear,  amplitude  is  arbitrary.  For  the  comparison,  the  amplitudes  of  the  PSE  and  DNS  density 
fluctuations  are  equated  at  the  inflow  station,  and  all  other  disturbance  quantities  are  scaled  propor¬ 
tionately.  Since  the  PSE  computation  is  performed  in  Fourier  space,  the  solution  provides  the 
contents  of  individual  modes  by  default.  To  extract  the  temporal  Fourier  harmonics  from  the  DNS 
data,  we  perform  a  time-series  analysis  analogous  to  hot-wire  anemometry  in  physical  experiments. 
The  flow  field  is  sampled  over  time  at  selected  grid  points  (usually  a  substantial  subset  of  the 
computational  grid)  and  is  subsequently  Fourier  analyzed  for  its  harmonic  content.  Figure  12  presents 
the  amplitude  envelope  of  the  density  component  of  the  fundamental  mode  versus  r],  obtained  from  32 
samples  that  span  the  one-period  interval  between  periods  42  and  43.  Because  the  phase  velocity  of 
the  second-mode  disturbance  is  about  90%  of  the  edge  velocity,  by  period  42  the  leading  wavefront 
has  exited  the  outflow  boundary  and  the  flow  has  settled  into  a  quasi-steady  periodic  state.  Individual 
profiles  in  Figure  12(a)  correspond  to  different  streamwise  locations,  equally  spaced  between  the 
inflow  boundary  and  a  point  32  wavelengths  downstream  (ahead  of  the  beginning  of  the  buffer 
region).  For  clarity,  these  profiles  are  each  staggered  by  a  factor  of  2  on  the  logarithmic  plot;  hence, 
their  maxima  are  relative,  not  absolute.  Figure  12(b)  graphs  the  density  maxima  obtained  from  Figure 
12(a)  versus  the  streamwise  coordinate  x.  In  Figure  13  the  density  maxima  and  the  maxima  of  other 
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flow  quantities  obtained  in  analogous  fashion  are  plotted  versus  Rej^  (rather  than  versus  x)  and  are 
compared  with  the  PSE  results.  For  clarity,  only  every  eighth  DNS  value  is  plotted.  As  shown  in 
Figure  13,  the  agreement  between  the  PSE  and  DNS  results  is  excellent.  The  rapid  divergence  of  the 
DNS  and  PSE  results  near  the  outflow  boundary  is  attributed  to  the  nonphysical  damping  of  the 
instability  wave  as  it  traverses  the  buffer  region  and  should  be  disregarded.  From  Figure  13,  it  can  be 
observed  that,  in  a  nonparallel  flow,  maxima  (or  minima)  of  the  various  components  of  the  distur¬ 
bance  do  not  necessarily  occur  at  the  same  streamwise  locations.  Hence,  neutral  points  and  growth 
rates  are  nonuniquely  defined.  Based  on  the  temperature  component  of  the  disturbance,  the  lower  and 
upper  branch  neutral  points  occur  at  approximately  Re^  =  740  and  Re^  =  1050,  respectively,  in 
surprisingly  good  agreement  with  Figure  10. 

Although  Figure  13  shows  good  agreement  at  maxima,  it  provides  no  information  about  the 
disturbance  structure.  In  Figure  14  we  compare  the  velocity  and  density  profiles  of  the  PSE  and  DNS 
calculations  at  the  station  Rci^  =  1046  near  the  location  of  the  global  maximum  of  the  temperature 
fluctuation.  For  comparison  of  the  disturbance  structure,  we  scale  the  DNS  and  PSE  results  so  that 
their  respective  density  maxima  are  unity.  Before  rescaling,  the  respective  maxima  differ  by  less  than 
0.5%.  For  clarity,  only  every  third  DNS  value  (denoted  by  symbols)  is  shown  in  the  figure.  Despite  the 
complex  nature  of  the  disturbance  structure,  the  agreement  is  excellent.  Moreover,  this  agreement 
confirms  that  the  rapid  change  in  the  disturbance  structure  at  the  three  stations  that  are  farthest 
downstream  in  Figure  12(a)  is  a  physical  change  and  not  the  result  of  a  computational  anomaly. 
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Figure  13.  Case  2:  PSE  and  DNS  disturbance  amplitude  max¬ 
ima  versus  Re,. 


Figure  14.  Case  2:  Disturbance  structure  at  =  1046. 
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Figure  15.  Case  3:  Amplitude  maxima  of  fundamental  of  tem¬ 
perature  fluctuation  and  first  two  harmonics. 
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Figure  16.  Case  3:  Structures  of  fundamental  and  first  two 
harmonics  of  temperature  fluctuation  at  =  1050. 


Case  3:  Flat  Plate,  Mach  4.5  Nonparallel  Flow,  Two-Dimensional  Nonlinear  Disturbance.  With  regard 
to  parameter  values,  this  case  is  identical  to  the  case  above,  except  that  the  forcing  amplitude  is  quite 
large  at  e  =  0.1.  At  this  level  of  forcing,  temperature  fluctuations  will  grow  in  amplitude  to  a  value 
that  is  25%  of  the  edge  temperature  of  the  base  flow.  For  the  DNS  calculation,  wall-normal 
resolution  remains  as  before;  however,  32  grid  points  per  primary  disturbance  wavelength  are  used  for 
a  total  of  1153  points  in  x.  Based  on  the  resolution  criterion  established  for  the  linear  case,  the 
fundamental  and  its  first  harmonic  are  well  resolved.  The  second  harmonic,  however,  is  somewhat 
underresolved.  Higher  harmonics,  which  contain  very  little  energy,  are  unresolved. 

•Figure  15  compares  the  PSE  and  DNS  results  with  regard  to  the  streamwise  evolution  of  the 
maxima  of  the  temperature  component  of  the  fundamental  and  its  first  two  harmonics.  For  the  DNS, 
the  Fourier  analysis  is  performed  on  32  time  samples  that  span  the  interval  between  periods  48  and 
49.  The  methods  agree  well  qualitatively  and  quantitatively;  however,  greater  disagreement  occurs 
as  the  index  of  the  harmonic  increases.  At  the  streamwise  location  of  maximum  disturbance  ampli¬ 
tude,  the  PSE  and  DNS  results  differ  by  0.15%,  6%,  and  12%,  respectively,  for  the  fundamental  (1,0), 
the  first  harmonic  (2, 0),  and  the  second  harmonic  (3, 0).  (Recall  that,  because  the  present  calculation  is 
two-dimensional,  =  0.) 

A  further  comparison  is  shown  in  Figure  16,  in  which  the  PSE  and  DNS  predictions  of  the 
structure  of  the  temperature  component  of  the  fundamental  and  its  first  two  harmonics  are  compared. 
In  contrast  to  the  (linear-amplitude)  results  presented  in  Figure  14,  the  present  (nonlinear-amplitude) 
results  have  not  been  rescaled  and  are  presented  on  a  logarithmic  scale  so  that  all  three  modes 
can  be  shown  with  clarity.  Between  the  wall  and  the  critical  layer,  both  methods  produce  similar 
amplitude  distributions.  Outside  the  boundary  layer,  there  is  divergence  between  the  PSE  and  DNS 
results  in  the  asymptotic  rate  of  decay  of  the  second  harmonic.  There  are  any  number  of  possible 
physical  and/or  numerical  reasons  for  this  discrepancy.  Among  these  are  the  differences  between 
PSE  and  DNS  in  the  imposition  of  far-field  boundary  conditions.  In  the  PSE  approach  Dirichlet 
conditions  are  imposed  on  the  disturbance  in  Fourier  space  at  very  large  t\  (typically  tj  »  50).  Details 
of  the  imposition  of  boundary  conditions  for  the  PSE  method  can  be  found  in  Chang  et  al.  (1991).  In 
the  DNS,  inviscid  characteristic  conditions  (described  previously)  are  imposed  in  physical  space  along 
the  far-field  boundary  r\  =  7.5.  However,  we  remind  the  reader  that  the  logarithmic  scale  of  Figure  16 
exaggerates  these  differences,  which  are  of  order  10“®  relative  to  the  base  state,  and  it  is  beyond  our 
current  intentions  to  ferret  out  the  sources  of  errors  of  this  magnitude  or  to  state  categorically  which 
method  is  “correct.”  Within  the  boundary  layer,  the  region  of  principal  interest,  the  methods  agree 
remarkably  well. 

Finally,  a  major  effect  of  the  strong  nonlinearity  is  a  mean  (0, 0)  component  that  distorts  the  base 
flow.  Figure  17  shows  that  the  PSE  and  DNS  calculations  both  predict  large  (5%)  distortions  of  the 
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Figure  17.  Case  3:  Nonlinearly  generated  mean  temperature 
distortion  versus  rj  at  =  1050. 


base  temperature  distribution  near  the  critical  layer.  Disagreement  between  the  methods  is  roughly 
6%  at  the  points  of  maximum  distortion.  The  reason  for  the  disparity  is  presently  unknown. 

Case  4:  Sharp  Cone,  Mach  6.8  Nonparallel  Flow,  Three-Dimensional  Linear  Disturbance.  In  this  final 
validation  case  we  consider  the  high-speed  boundary-layer  flow  along  an  axisymmetric  sharp  cone, 
which  is  forced  to  excite  two  symmetric  oblique  (helical)  second-mode  disturbances.  The  geometry  and 
the  flow  parameters  given  below  are  chosen  to  correspond  approximately  to  the  stability  experiment 
of  Stetson  et  al.  (1983)  in  a  hypersonic  wind  tunnel  with  a  free-stream  Mach  number  of  8: 

cp  =  7°  (constant),  =  6.8,  T*  =  71  K.  (36) 

We  assume  that  the  boundary-layer  edge  conditions  remain  constant,  a  reasonable  approximation 
except  near  the  tip  of  the  cone,  a  region  excluded  in  the  present  simulations.  The  PSE  calculation 
spans  from  700  <  ^  1800.  To  limit  the  size  of  the  DNS  computation,  we  consider  only  the  region 

corresponding  to  1013  <  <  1536.  At  the  DNS  inflow  boundary,  the  disturbance  parameters 

are 

F  =  2.094  X  10-^ 

n  =  13  ipo  =  1.1736,  Ro  =  11.077),  (37) 

e  =  0.001, 

and  S^JL*  =  11.17.  Whereas  the  geometry  and  flow  parameters  correspond  to  the  experiment  of 
Stetson  et  al.  (1983),  the  disturbance  parameters  were  selected  somewhat  arbitrarily.  The  dimensionless 
frequency  given  in  (37)  corresponds  to  a  physical  frequency  of  180  kHz,  considerably  higher  than  the 
dominant  102  kHz  frequency  observed  in  the  experiment.  The  higher  frequency  was  selected  for  this 
validation  case  because  its  entire  instability  region,  including  the  lower  and  upper  branch  neutral 
points,  is  contained  within  a  fairly  short  stream  wise  region,  one  of  (almost)  reasonable  size  for  a  DNS 
calculation  of  modest  computational  effort.  There  is  nothing  sacred  about  the  particular  choice  of 
n  =  13,  except  that  it  yields  an  unstable  mode  of  a  sizeable  obliqueness  angle  (24.66°  at  the  DNS 
inflow  boundary)  for  the  given  frequency. 

Figure  18  presents  the  evolution  of  the  amplitudes  of  each  of  the  five  components  of  the  fundamen¬ 
tal  (1, 1)  mode  for  700  <  Re^  <  1500.  In  this  case  the  buffer  region  for  the  DNS  calculation  is  not 
shown.  As  in  Case  2,  the  PSE  calculation  is  performed  in  the  linear  mode,  in  which  case  amplitude  is 
arbitrary.  For  the  comparisons  below,  the  maximum  densities  of  the  PSE  and  DNS  calculations 
are  equated  at  the  DNS  inflow  plane,  and  all  other  quantities  are  scaled  accordingly.  The  agreement 
between  the  DNS  and  PSE  results  is,  in  general,  quite  good.  The  kinks  in  the  u'  maxima  curve  at 
=  1120  and  the  w'  maxima  curve  at  Re^  =  1280  are  due  to  the  double-humped  structure  of  the 
disturbance  shown  in  Figure  19.  Remarkably,  the  DNS  and  PSE  results  predict  the  shift  of  the  global 
maximum  from  one  hump  to  the  other  at  precisely  the  same  values  of  Re^.  There  is  some  discrepancy 
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Figure  18.  Case  4;  PSE  and  DNS  disturbance  amplitude  maxima  versus  Re^  for  cone  case;  (a)  T  and  p'  and  (b)  u',  v',  and  w'. 
Two  distinct  PSE  results  differ  only  in  the  method  by  which  wall-normal  derivatives  of  the  base  flow  are  computed. 


in  the  DNS  and  PSE  maxima,  which  is  most  significant  in  the  temperature,  and  which  increases  with 
stream  wise  distance.  To  eliminate  one  potential  source  of  disagreement  for  this  comparison,  the 
laminar  base  flow  Q°  along  the  cone  was  obtained  for  both  the  PSE  and  DNS  calculations  from  the 
spectrally  accurate  BE  code  of  Pruett  and  Streett  (1991)  and  Pruett  (1993).  The  remaining  sources  of 
disagreement  may  be  either  physical  or  numerical  or  a  combination.  We  remind  the  reader  that  the 
governing  equations,  the  far-field  boundary  conditions,  and  the  discretions  of  the  two  methods  differ, 
to  mention  but  a  few  of  many  subtle  differences  whose  effects  may  be  cumulative.  Because  the  flow  is 
sensitive  to  seemingly  imperceptible  changes,  it  is  somewhat  pointless  to  try  to  isolate  the  sources  of 
disagreement.  Perhaps  it  is  more  instructive  to  demonstrate  this  flow-field  sensitivity.  In  Figure  18  two 
sets  of  PSE  results  are  shown.  The  two  sets  differ  only  in  the  method  by  which  derivatives  of  the  base 
state  were  obtained.  In  the  first  set  all  derivatives  were  computed  internally  in  the  PSE  code  to 
fourth-order  accuracy.  In  the  second  set  the  derivatives  u^,  T^,  and  were  computed  to 

spectral  accuracy  in  the  BE  code  and  were  subsequently  passed  to  the  PSE  code.  Figure  20  compares 
the  first  and  second  derivatives  of  the  temperature  as  computed  by  the  BE  and  PSE  codes.  At  worst, 
the  temperature  derivatives  computed  by  the  two  methods  differ  by  less  than  0.001,  and,  in  the 
graphical  format,  they  are  indistinguishable.  Moreover,  velocity  derivatives  obtained  by  the  two 
methods  agree  everywhere  to  more  than  four  significant  digits.  Nevertheless,  the  cumulative  effect  of 
these  seemingly  insignificant  differences  between  derivative  profiles  is  discernible  in  Figure  18.  This 
inconsistency  did  not  arise  for  the  flat-plate  test  cases  studied  previously,  for  which  the  base  flow  is 


Figure  19.  Case  4:  Structure  at  Re^  =  1291  of  (a)  T'  and  p'  and  (b)  u',  v',  and  w'. 


74 


C.D.  Pruett,  T.A.  Zang,  Chau-Lyan  Chang,  and  M.H.  Carpenter 


0.5 


0.0 


-0.5 


-1.0 


0.0  0.5  1.0  1.5  2.0 

Figure  20.  Case  4:  Wall-normal  derivatives  of  base  tempera- 
7y  ture. 


obtained  from  a  similarity  solution  whose  derivatives  can  be  obtained  semianalytically.  These  results 
reveal  the  extreme  sensitivity  of  the  stability  of  high-speed  wall-bounded  flows  to  changes  in  the  base 
state  and  underscore  the  necessity  of  numerical  methods  of  the  highest  accuracy  for  such  problems. 

Before  concluding,  we  offer  a  few  comments  on  the  computational  resources  required  by  the 
present  algorithm.  DNS  is,  in  general,  an  expensive  research  tool.  Cases  3  and  4  of  this  section 
consumed  approximately  40  and  150  hours,  respectively,  on  a  single  processor  of  a  Cay  Y-MP.  It  is 
estimated  that  spatial  DNS  of  the  complete  laminar-breakdown  process  for  a  high-speed  boundary- 
layer  flow  will  require  at  least  2000  Cray-2  hours  and  at  least  256  megawords  of  memory.  Needless  to 
say,  such  a  tool  is  not  for  routine  use  at  the  present  time.  Nevertheless,  spatial  DNS  is  invaluable  as  a 
means  of  validating  less-expensive  approximate  methods  (e.g.,  PSE  and  LES)  and  for  building 
high-fidelity  data  bases  of  transitioning  flows,  from  which  transition  modelers  can  construct  simplified 
models.  With  continued  algorithm  refinements,  supercomputer  advancements,  and  massively  parallel 
impelementations,  we  believe  spatial  DNS  of  transition  to  turbulence  will  be  practical  and  routine 
well  within  a  decade. 


6.  Conclusions 

A  highly  accurate  algorithm  has  been  developed  for  the  direct  numerical  simulation  (DNS)  of  forced, 
spatially  evolving  instability  waves  in  high-speed  wall-bounded  flows.  To  minimize  dissipation  and 
dispersion  errors,  the  fully  explicit  algorithm  exploits  both  spectral  collocation  and  high-order  central 
compact-difference  techniques.  In  its  present  form  the  algorithm  allows  for  three-dimensional  flow 
along  two-dimensional  or  axisymmetric  bodies.  Of  particular  interest  in  this  work  are  the  flat-plate 
and  the  sharp-cone  geometries. 

Part  1  of  this  paper  deals  primarily  with  thorough  validation  of  the  DNS  scheme  by  comparisons 
with  results  obtained  from  classical  linear  stability  theory  (LST)  and  from  the  parabolized  stability 
equation  (PSE)  method.  Test  cases  examine  forced  two-dimensional  second-mode  instability  waves  in 
Mach  4.5  flat-plate  boundary-layer  flows  and  three-dimensional  second-mode  waves  in  a  Mach  6.8 
flow  along  a  sharp  cone.  From  these  validation  studies,  several  insights  emerge  that  pertain  both  to 
the  numerical  methods  and  to  the  physical  problems  addressed.  First,  for  the  Mach  numbers  of  the 
test  cases,  the  streamwise  resolution  needed  to  obtain  good  agreement  with  theoretical  results 
was  higher  than  anticipated.  With  sixth-order  streamwise  accuracy,  approximately  12  and  16  grid 
points  per  wavelength  were  required  in  order  to  extract  the  growth  rate  of  a  monochromatic 
disturbance  to  three  and  four  place  accuracy,  respectively.  Second,  for  nonparallel  flows,  consistent 
inflow  conditions  for  the  DNS  can  be  obtained  from  the  PSE  method,  whereas  inflow  conditions 
derived  from  LST  are  fundamentally  inconsistent  due  to  the  parallel-flow  approximation  of  LST. 
Third,  at  these  high  Mach  numbers,  the  stability  of  the  boundary  layers  was  found  to  be  sensitive  to 
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the  slightest  (almost  imperceptible)  change  in  the  base  state.  Fourth,  for  low-amplitude  (linear) 
disturbances  the  agreement  between  the  PSE  method  and  the  DNS  results  was  near  perfect  in  terms 
of  both  the  amplitude  and  the  structure  of  the  disturbances.  Fifth,  for  large-amplitude  (nonlinear) 
disturbances,  agreement  between  the  two  methods  was  near-perfect  for  the  fundamental  and  good  for 
engineering  purposes  for  the  higher  harmonics,  with  the  tendency  towards  greater  disagreement  the 
higher  the  harmonic.  Finally,  the  generally  close  agreement  between  DNS  and  PSE  results  convinc¬ 
ingly  demonstrates  the  potential  of  the  PSE  method  as  a  reliable  new  tool  for  analysis  of  hydro- 
dynamic  stability  in  high-speed  wall-bounded  flows. 

At  present  there  is  a  great  need  for  experimental  studies  of  stability  and  transition  in  high-speed 
boundary-layer  flows.  Until  such  experiments  are  accomplished,  the  only  available  source  of  detailed 
flow-field  information  for  high-speed  transitional  flows  is  DNS.  The  development  and  validation  of 
the  present  DNS  algorithm  thus  fills  a  void  in  the  national  capability  in  the  area  of  transition 
research.  Using  this  new  tool,  Part  2  of  this  work  will  examine  in  detail  the  laminar  breakdown  of  a 
perturbed  high-speed  boundary-layer  flow  along  a  cone. 
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Abstract 

A  suite  of  highly  accurate  numerical  algorithms  has  been  developed  and  validated 
for  use  in  the  simulation  of  crossflow  instabilities  on  supersonic  swept  wings,  an  appli¬ 
cation  of  potential  relevance  to  the  design  of  the  High-Speed  Civil  Transport  (HSCT). 
Principal  among  these  algorithms,  and  the  primary  focus  of  this  report,  is  a  direct 
numerical  simulation  (DNS)  scheme  embodied  in  the  code  CMPSBL,  which  solves 
the  unsteady,  three-dimensional  compressible  Navier-Stokes  equations  in  orthogonal, 
body-fitted  coordinates.  The  DNS  algorithm  is  fully  explicit  in  time  and  exploits  a 
combination  of  spectral  and  high-order  compact- difference  techniques  for  spatial  dis¬ 
cretizations.  A  companion  code  WINGBL2,  documented  in  a  previous  technical  report, 
exploits  spectral- collocation  techniques  to  solve  the  compressible  boundary-layer  equa¬ 
tions  to  provide  an  accurate  basic  state  to  the  DNS.  In  addition,  the  algorithm  INTBL 
uses  spectral  techniques  to  interpolate  the  boundary-layer  solution  onto  the  computa¬ 
tional  grid  of  the  DNS.  These  algorithms  are  then  used  to  examine  the  development  of 
stationary  crossflow  instability  on  an  infinitely  long  77-degree  swept  wing  in  Mach  .3.5 
flow.  Crossflow  disturbances  are  generated  by  simulating  spanwise-periodic  roughness 
elements  downstream  of  the  computational  inflow  boundary.  The  results  of  the  DNS  are 
compared  with  the  predictions  of  linear  parabolized  stability  equation  (PSE)  method¬ 
ology  obtained  with  the  code  ECLIPSE  (developed  independently  of  this  effort).  In 
general,  the  DNS  and  PSE  results  agree  closely,  thereby  providing  a  reasonable  valida¬ 
tion  of  both  approaches.  Specifically,  both  methods  show  the  alignment  of  stationary 
crossflow  vortices  along  inviscid  streamlines  as  anticipated.  Moreover,  the  methods 
agree  well  in  the  predicted  spatial  evolution  of  a  small-amplitude  (linear)  crossflow 
mode  and  in  the  structure  of  this  mode.  Although  further  validation  is  warranted  (for 
large- amplitude  stationary  and  traveling  crossflow  disturbances),  the  present  results 
demonstrate  a  new  numerical  capability  relevant  to  a  problem  of  practical  importance. 
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1  Introduction 


Understanding,  predicting,  and  controlling  laminar-turbulent  transition  remains  the  holy 
grail  of  aeronautics  research  despite  more  than  a  century  of  assault  from  the  combined  forces 
of  theory,  experiment,  and  computation.  Recent  advances  on  each  of  these  fronts  and  in  the 
areas  of  materials  processing  and  micro- actuators  (Ho  and  Tai  [13]),  however,  have  brought 
this  elusive  goal  within  sight,  and  there  is  renewed  interest  in  laminar-flow  control  (LFC) 
technology. 

Because  surface  friction  and  heat  transfer  increase  dramatically  as  the  boundary  layer 
transitions  from  a  laminar  to  a  turbulent  state,  the  design  of  efficient  aerospace  vehicles 
depends  upon  accurately  predicting  the  transition  location  and  the  extents  of  the  regions  of 
laminar,  transitional,  and  turbulent  flows.  Although,  in  the  current  economic  environment, 
it  is  questionable  whether  laminar  flow  control  (LFC)  technology  can  be  made  commercially 
viable  in  the  near  future  for  subsonic  aircraft,  the  potential  economic  benefit  could  be  sig¬ 
nificant  for  supersonic  transport  aircraft  such  as  the  proposed  High-Speed  Civil  Transport 
(HSCT).  This  is  because  maintenance  of  laminar  flow  over  a  substantial  portion  of  the  wing 
of  the  HSCT,  for  example,  would  not  only  reduce  drag  but  would  also  reduce  thermal  loads 
on  the  structure. 

The  practical  attainment  of  LFC  technology  will  require  fundamental  understanding 
of  the  stability  of  three-dimensional  boundary-layer  flows.  Examples  of  prototypical  three- 
dimensional  boundary-layer  flows  are  flows  past  rotating  cones  and  spheres,  flow  on  a  rotating 
disk,  flows  in  corners,  and  flows  over  sw’ept  wings,  the  latter  of  which  is  of  the  most  practical 
relevance. 

Whereas  the  stability  of  two-dimensional  flows  has  been  studied  extensively,  researchers 
are  only  beginning  to  focus  attention  on  three-dimensional  boundary- layer  flows.  Theoreti¬ 
cal  and  experimental  results  demonstrate  that  there  exists  a  far  greater  variety  of  paths  to 
transition  for  three-dimensional  boundary  layers  than  for  two-dimensional  flows.  For  exam¬ 
ple,  in  the  case  of  the  swept  wing,  the  flow  may  be  susceptible  to  Tollmien-Schlicting  (TS) 
instability,  Goertler  instability  (associated  with  the  concavity  of  the  lower  surface;  e.g.,  see 
Hall  [12]),  crossflow  instability  (the  subject  of  the  discussion  to  follow),  and  attachment-line 
instability  (instability  of  the  flow  along  the  leading  edge;  e.g.,  see  Joslin  [17]).  For  the  wing 
of  the  HSCT,  to  be  specific,  crossflow  instability  is  likely  to  be  the  most  "dangerous’'  in  the 
sense  of  leading  most  rapidly  to  transition. 

Crossflow  instability  w’as  first  identified  in  the  1940’s  in  experiments  related  to  the 
Northrop  flying  wing.  As  stated  by  Reed  and  Saric  [28]  in  their  review  paper,  “Although 
unyawed  wind-tunnel  tests  showed  laminar  flow  back  to  60  percent  chord,  yawed  flight  tests 
showed  turbulent  flow  from  the  leading  edge  on  both  the  upper  and  lower  surfaces.  ’  In  flow 
over  a  swept  wing,  the  inviscid  streamlines  form  S-curves  because  the  pressure  gradients 
associated  with  body  curvature  accelerate  or  decelerate  the  flow  in  the  direction  perpendic¬ 
ular  to  the  leading  edge  while  leaving  the  velocity  component  parallel  to  the  leading  edge 
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virtually  unchanged.  The  combination  of  curvilinear  inviscid  streamlines  and  the  viscous 
no-slip  condition  at  the  wall  generates  a  crossflow  velocity  perpendicular  to  the  local  inviscid 
streamline.  The  crossflow  velocity  is  zero  both  at  the  wall  and  in  the  freestream  (by  defini¬ 
tion),  and  thus,  in  between  it  experiences  a  maximum  and  a  point  of  inflection.  Maxima  for 
crossflow  velocities  are  typically  on  the  order  of  3-4  percent  of  the  velocity  in  the  direction 
of  the  streamlines. 

Even  though  crossflow  velocities  are  typically  relatively  small,  the  presence  of  the  in¬ 
flection  point  renders  the  flow  susceptible  to  crossflow  instability,  which  is  of  inviscid  type. 
The  instability  manifests  itself  as  co-rotating  vortices  that  align  themselves  roughly  along 
inviscid  streamlines.  The  spacing  between  adjacent  vortices  is  on  the  order  of  several  S*,  the 
boundary- layer  displacement  thickness.  According  to  Mack  [20],  the  instability  exists  for 
a  whole  band  of  frequencies,  including  zero.  Curiously,  linear  stability  theory  predicts  the 
most  amplified  disturbances  to  be  traveling  waves,  whereas  stationary  crossflow  instability 
is  frequently  observed  in  experiments,  except  close  to  laminar  breakdown  (Reed  and  Saric 
[28]  ).  Choudhari  and  Streett  [10]  and  Choudhari  [9]  have  shown,  on  the  basis  of  receptivity 
theory,  that  surface  irregularities  may  favor  the  development  of  stationary  crossflow  modes 
by  giving  them  much  larger  initial  amplitudes  than  those  of  traveling  crossflow  waves. 

Both  theory  and  experiment  concur  that  crossflow  instability  dominates  on  swept  wings  in 
regions  of  rapidly  changing  pressure.  Other  regions  on  the  wing,  however,  may  be  dominated 
by  TS-like  instabilities  (which  we  take  to  include  the  first-mode  instability  of  compressible 
boundary-layer  flow).  According  to  Reed  and  Saric  [28],  a  major  unanswered  question  is  "the 
interaction  between  crossflow  vortices  and  TS  waves.”  It  appears  that  crossflow  vortices  can 
modify  the  TS  instability  to  enhance  its  growth  rate.  From  experiments,  it  also  appears 
that  the  crossflow  instability  is  extremely  sensitive  to  initial  conditions,  and  it  has  even  been 
suggested  that  the  theory  for  crossflow  instability  is  not  well  posed  (Reed  and  Saric  [28]). 

In  summary,  transition  on  a  swept  wing  is  a  highly  sensitive  and  complex  process.  The 
transition  location  is  affected  by  nonparallelism  of  the  mean  flow,  pressure  gradient,  surface 
roughness,  freestream  turbulence,  and  body  curvature.  Moreover,  these  elements  give  rise  to 
both  inviscid  crossflow  modes  and  to  TS-like  instabilities,  which  may  interact.  It  is  unrea¬ 
sonable  to  ask  any  approximate  theory  to  accommodate  all  of  these  diverse  and  interacting 
elements.  As  a  result,  the  problem  is  well  suited  for  numerical  investigations  using  direct 
numerical  simulation  (DNS),  for  which  a  minimum  of  simplifying  assumptions  are  made. 

Crossflow  instability  on  swept  wings  in  incompressible  flows  has  been  investigated  re¬ 
cently  with  a  parabolized  stability  equation  (PSE)  approach  by  Malik  et  al.  [21]  and  with 
DNS  by  Lin  and  Reed  [19],  Fuciarelli  and  Reed  [11],  Joslin  and  Streett  [18],  and  Joslin 
[16].  The  present  work  is  believed  to  be  the  first  investigation  of  crossflow  instability  on  a 
supersonic  swept  wing  by  means  of  DNS.  Direct  simulation  is  computationally  expensive; 
consequently,  we  believe  that  simulations  are  most  fruitful  when  focused  along  with  theoret¬ 
ical  and  experimental  investigations  on  a  single  problem  of  practical  importance.  To  that 
end,  we  have  selected  a  test  problem  that  parallels  the  quiet  wind-tunnel  experiment  on  a 
supersonic  swept  wing  by  Cattafesta  et  al.  [3].  Moreover,  our  DNS  results  are  compared 


with  results  obtained  by  PSE  methodology  for  compressible  flows  (e.g.,  Chang  et  al.  [6]) 
in  an  attempt  to  cross-validate  both  methods  for  this  difficult  problem.  For  computational 
efficiency,  we  must  limit  consideration  to  quasi-three-dimensional  (infinite-span)  swept-wing 
flows.  This  assumption  permits  highly  efficient  spectral  collocation  methods  to  be  exploited 
in  the  spanwise  direction.  The  limitation  of  the  computational  model  should  be  kept  in 
mind  in  drawing  conclusions  relative  to  the  experiment.  We  emphasize,  however,  that  the 
full  effects  of  streamwise  surface  curvature  are  incorporated  into  the  DNS  model,  in  contrast 
to  most  recent  approaches. 

The  next  section  discusses  the  coordinate  system  and  governing  equations.  Section  3 
summarizes  the  numerical  methodology.  The  computational  test  case  is  addressed  in  Section 
4.  Results  are  presented  in  Section  5,  and  Section  6  concludes  with  a  few  closing  remarks. 


2  Governing  Equations 


The  flow  is  governed  by  the  compressible  Navier-Stokes  equations  in  the  form  given  in  Eqs. 
(1)-(15)  of  Pruett  et  al.  [26],  which  is  appropriate  for  body-fitted  coordinates  on  either 
axisymmetric  or  two-dimensional  bodies.  Here,  we  consider  only  two-dimensional  (infinite 
in  span)  wing-like  bodies.  Let  the  wing  be  imbedded  in  a  rectangular  cartesian  coordinate 
system  (i^,  r},  C)?  with  the  wing  cross  section  specified  as  C(^).  Let  (x,  y,  z)  denote  body-fitted 
coordinates  such  that  x  is  the  surface  arc  length  normal  to  the  attachment  line,  y  is  the 
spanwise  coordinate  (parallel  to  the  attachment  line),  and  2  is  the  wall-normal  coordinate. 
(The  reader  will  note  that,  for  consistency  with  the  coordinates  of  Pruett  et  al.  [26],  the 
y  and  ^  coordinates  are  switched  relative  to  the  normal  convention.)  Accordingly,  u,  u, 
and  10  are  the  velocity  components  in  the  streamwise,  spanwise,  and  wall-normal  directions, 
respectively.  The  fundamental  metric  tensor  has  only  one  non-constant  quantity  s,  which 
arises  from  streamwise  curvature  and  is  defined  in  Eq.  (2)  of  Pruett  et  al.  [26].  Moreover, 
by  p,  p,  r,  p,  and  k,  we  denote  the  density,  pressure,  temperature,  viscosity,  and  thermal 
conductivity  of  the  fluid,  respectively.  The  viscosity  p  is  modeled  by  Sutherland’s  law, 
and  K  =  p/Pr,  where  Pr  is  the  Prandtl  number.  For  computational  efficiency,  the  energy 
equation  is  cast  in  terms  of  the  pressure,  as  given  in  Eq.  (15)  of  Pruett  et  al.  [26].  All  flow 
variables,  except  pressure,  are  non-dimensionalized  by  post-shock  reference  values,  denoted 
by  subscript  r.  Pressure  is  scaled  by  Lengths  are  scaled  by  the  boundary-layer 

displacement  thickness  8*  at  the  (computational)  inflow  boundary.  Throughout  this  work, 
dimensional  quantities  are  denoted  by  asterisk. 


3  Methodology 


The  conventional  approach  to  (spatial)  stability  analysis,  adopted  here,  consists  of  three 
basic  steps:  determination  of  a  laminar  base  state  whose  stability  is  to  be  investigated. 
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perturbation  of  the  base  state  by  superposition  of  disturbances  at  or  near  the  computational 
inflow  boundary,  and  calculation  of  the  spatial  evolution  of  the  disturbances.  Each  of  these 
steps  is  addressed  in  turn  below. 


3.1  Base  Flow 

The  base  state  is  computed  by  the  spectrally  accurate  boundary- layer  code  VVINGBL2  of 
Pruett  [25],  which  was  specifically  designed  for  the  infinite-span  swept-wing  problem.  The 
effects  of  streamwise  curvature,  streamwise  pressure  gradient,  and  wall  suction/blowing  are 
taken  into  account  in  the  governing  equations  and  boundary  conditions.  The  boundary-layer 
equations  are  formulated  both  for  the  attachment-line  flow  and  for  the  evolving  boundary 
layer.  The  equations  for  the  evolving  flow  are  solved  by  an  implicit  marching  procedure 
in  the  direction  perpendicular  to  the  leading  edge,  for  which  high-order  (up  to  5th)  back¬ 
ward  differencing  techniques  are  used.  In  the  wall-normal  direction,  a  spectral  collocation 
method  based  on  Chebyshev  polynomial  approximations  is  exploited.  Spectral  accuracy  is 
advantageous  in  that  1)  the  solution  is  highly  accurate  even  for  relatively  coarse  grids,  2)  the 
boundary-layer  profiles  and  their  derivatives  are  extremely  smooth  (a  necessity  for  stability 
analyses),  and  3)  interpolation  to  other  grids  can  be  accomplished  with  virtually  no  loss  of 
accuracy. 

A  few  comments  in  regard  to  point  3)  above  are  in  order.  Because  WINGBL2  is  de¬ 
signed  specifically  for  applications  to  stability  analyses,  DNS,  and  large-eddy  simulation 
(LES),  special  attention  has  been  paid  to  the  process  of  interpolating  data  to  other  grids, 
for  example,  to  a  DNS  grid.  For  this  purpose,  a  companion  code  INTBL  was  written.  In 
brief  the  procedure  is  as  follows.  In  the  output  from  WINGBL2,  lengths  are  scaled  by  d*. 
which  grows  with  x*.  The  outer  edge  of  the  boundary  layer  typically  lies  between  2  and  4 
displacement  thicknesses  from  the  wall.  Let  5re(a^)  denote  the  location  of  the  boundary-layer 
edge.  To  interpolate  to  a  grid  uniform  in  2,  for  example,  we  first  perform  spectrally  accurate 
interpolation  in  the  interior  region  0  <  2  <  re  coupled  with  analytic  extrapolation  outside 
the  boundary  layer  (z  >  Ze).  This  step  is  followed  by  high-order  (typically  5th)  polynomial 
interpolation  in  x. 

The  analytic  extrapolation  in  the  far  field  is  accomplished  at  each  x  by  solving  the 
ordinary  differential  equation  that  results  from  the  continuity  equation  in  the  asymptotic 
limit  2  ^  00.  In  many  boundary-layer  approximations,  some  curvature  effects  are  neglected, 
in  which  case  the  continuity  equation  may  be  inexact.  Indeed,  it  was  determined  that  Eq. 
(7-100)  on  page  397  of  Anderson  et  al.  [1],  on  which  the  continuity  equation  for  WINGBL2 
was  originally  based,  is  inexact  for  compressible  flow;  an  exact  continuity  equation  was 
subsequently  derived.  That  the  continuity  equation  be  exact  is  essential  if  the  near-field 
solution  and  the  far-field  extrapolation  are  to  merge  continuously,  as  they  must  for  DNS. 
The  method  proposed  by  Pruett  [24]  to  extract  wall-normal  velocity  accurately,  adopted  here 
also,  permits  an  independent  check  of  continuity.  Typically  WINGBL2  conserves  continuity 
nearly  to  machine  precision,  as  shown  in  Figs.  1  and  2. 
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The  input  to  WINGBL2  consists  primarily  of  three  files  that  contain  the  wing  geometry, 
the  wing  pressure  distribution  in  the  form  of  tabulated  pressure  coefficients,  and  the  reference 
conditions.  Cubic  splines  are  used  to  interpolate  as  necessary  between  the  tabulated  values. 
The  reference  conditions  are  presumed  to  be  downstream  of  leading-edge  shocks. 

3.2  Disturbances 

The  base  flow  is  perturbed  in  the  manner  of  Joslin  [16],  who  simulated  crossflow  instability 
in  incompressible  flow  on  a  swept  wing.  Disturbances  are  introduced  by  mass-preserving 
suction  and  blowing  at  the  wall.  Currently,  we  exploit  steady  suction/blowing  in  order 
to  induce  stationary  crossflow  modes.  The  suction/blowing  strip  spans  the  width  of  the 
wing,  but  is  localized  in  x.  (Henceforth,  for  brevity,  we  will  refer  only  to  the  “suction" 
strip.)  As  noted  by  Joslin  [16]:  “This  mode  of  disturbance  generation  would  correspond 
to  an  isolated  roughness  element  within  the  computational  domain.”  Joslin  has  found  by 
numerical  experimentation  that  the  suction  strip  should  not  be  located  farther  upstream  than 
approximately  10  percent  chord.  Otherwise,  the  crossflow  modes  are  not  swept  downstream. 
The  spanwise  wavenumber  jS  of  the  disturbance  strip  is  a  fundamental  parameter  of  the  flow, 
as  is  the  maximum  amplitude  of  the  wall  velocity  u^walh  Typically,  very  small  normalized 
wall  velocities  (say,  ir^ali  =  10“'*  or  less)  suffice  to  trigger  stationary  crossflow  instability. 
The  suction/blowing  profile  is  a  full  sine  wave  in  the  spanwise  direction;  to  ensure  continuity 
of  derivatives,  the  cube  of  a  half  sine  function  is  used  to  shape  the  wall  velocity  in  the 
streamwise  direction. 


3.3  DNS  Methodology 


For  the  present  application,  we  exploit  the  highly  accurate  DNS  methodology  of  Pruett  et 
al.  [26]  with  minor  modifications.  To  summarize  briefly,  time  is  advanced  fully  explicitly  by 
the  third-order  low-storage  Runge-Kutta  scheme  of  Williamson  [30].  Spatial  derivatives  are 
approximated  by  a  combination  of  spectral-collocation  techniques  and  high-order  compact- 
difference  schemes.  In  the  streamwise  and  wall-normal  directions,  we  exploit  fourth-  and 
sixth-order  compact-difference  operators,  respectively.  In  the  spanwise  direction,  the  flow 
is  assumed  to  be  periodic,  making  possible  the  use  of  a  spectral  collocation  technique  with 
Fourier  exponential  basis  functions.  The  assumption  of  spanwise  periodicity  restricts  the 
basic  flows  under  consideration  to  quasi-three-dimensional;  that  is,  all  base-flow  quantities, 
including  the  three  velocity  components,  are  assumed  to  be  functions  only  of  the  streamwise 
and  wall-normal  coordinates.  This  is  equivalent  to  assuming  that  the  wing  is  of  infinite  span 
with  constant  cross-section. 

The  boundary  conditions  require  some  modification  from  those  given  in  Pruett  et  al.  [26]. 
For  a  well-designed  wing,  the  streamwise  component  of  velocity  should  be  everywhere  sub¬ 
sonic.  Therefore,  the  streamwise  velocity  at  the  computational  inflow  boundary  is  subsonic. 
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and  in  the  outer  (Euler)  region  of  the  domain,  there  exists  an  upstream  characteristic  veloc¬ 
ity.  Accordingly,  we  specify  the  Riemann  invariants  along  the  inflow  boundary.  This  inflow 
treatment,  correct  for  inviscid  flow,  is  not  entirely  satisfactory  in  the  viscous  layer,  where  all 
flow  quantities  should  be  specified  (Poinsot  and  Lele  [22]).  The  wall  and  far-field  boundary 
treatments  are  the  same  as  described  in  Pruett  et  al.  [26],  except  that  suction  and  blowing 
are  introduced  at  the  wall  to  induce  the  disturbance,  as  described  previously.  We  currently 
exploit  a  buffer  domain  (Streett  and  Macaraeg  [29])  in  the  vicinity  of  the  outflow  boundary, 
as  was  done  successfully  by  Pruett  et  al.  [26].  However,  for  this  particular  application,  we 
are  presently  experiencing  some  reflection  from  the  outflow  boundary.  It  is  not  yet  known 
whether  the  reflection  is  of  numerical  or  physical  origin  (or  both).  Additional  effort  should 
be  directed  to  diagnosis  and  refinement  of  the  inflow  and  outflow  boundary  conditions. 

The  work  of  Pruett  et  al.  [26]  and  the  present  work  differ  conceptually  in  that,  for  the 
latter,  which  is  concerned  with  stationary  crossflow  instability,  only  the  time-asymptotic 
solution  is  of  interest.  Time  integration  is  accomplished  solely  as  a  means  of  relaxation 
toward  the  steady  state.  Physically,  from  the  point  at  which  a  disturbance  is  introduced,  a 
wav'e  packet  propagates  (predominantly)  downstream,  depositing  in  its  wake  the  stationary 
crossflow  instability  of  interest.  Once  the  leading  wave  packet  has  exited  the  domain,  the  flow- 
settles  to  its  perturbed  steady  state.  Attainment  of  a  steady  state  is  assessed  by  computing 
the  residuals  of  the  time-independent  compressible  Navier-Stokes  equations.  In  practice, 
in  the  context  of  fully  explicit  time  advancement,  the  residual  is  simply  the  discrete  update 
vector.  Evolution  of  the  global  maximum  residual  for  the  calculation  of  the  Results  section 
is  shown  in  Fig.  3.  The  largest  residuals  are  associated  with  the  continuity  [p]  and  spanwise 
momentum  (u)  equations.  The  maximum  residuals  grow  approximately  exponentially  in  time 
as  the  leading  wavefront  propagates  downstream  at  a  nearly  constant  velocity,  as  implied  by 
Fig.  4.  Different  formulations  for  the  buffer  domain  may  change  the  velocity  of  propagation 
in  this  region.  For  example,  the  buffer  domain  treatment  that  resulted  in  the  least  reflection 
also  unfortunately  significantly  diminished  the  propagation  velocity  in  the  buffer  region, 
necessitating  a  relatively  long  integration  time.  Following  the  exit  of  the  leading  wavefront 
from  the  domain  (not  shown),  the  residuals  decay  rapidly,  provided  the  outflow  boundary 
treatment  is  non-reflecting.  Moreover,  examination  of  the  local  residual  field  (also  not  shown) 
indicates  that  the  flow  is  essentially  stationary  a  short  distance  upstream  of  the  trailing  edge 
of  the  propagating  wave  packet.  Finally,  Fig.  5  shows  the  location  of  the  maximum  residual 
in  terms  of  the  wall-normal  gridpoint  index  k.  The  initial  outward  movement  of  the  residual 
is  associated  with  a  weak  acoustic  pulse  that  radiates  from  the  receptivity  region  near  the 
suction  strip.  Following  the  exit  of  this  wave  from  the  upper  boundary,  the  maximum 
residual  remains  near  k  —  40,  a  distance  of  approximately  1.5<5  from  the  wall,  until  the 
wave  encounters  the  buffer  domain.  The  wall-normal  location  of  maximum  residual  coincides 
closely  with  the  maximum  perturbation  amplitude  of  the  crossflow-  instability. 
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3.4 


PSE  Methodology 


The  PSE  method,  developed  originally  for  incompressible  flow,  has  been  extended  to  com¬ 
pressible  boundary  layers  by  Bertolotti  and  Herbert  [2]  and  by  Chang  and  coworkers  (See 
references  [4]  and  [7]).  The  method  is  rapidly  gaining  favor  as  a  powerful  and  efficient  tool 
for  analyzing  the  stability  of  spatially  evolving  boundary  layers.  The  method  treats  both 
nonparallel  boundary-layer  effects  and  moderately  nonlinear  wave  interactions.  The  method 
has  recently  been  adapted  for  application  to  supersonic  swept  wings  by  Chang  and  coworkers. 
The  theory  is  presented  in  Reference  [6];  a  user-friendly  code  ECLIPSE  for  industry  applica¬ 
tions  has  also  been  developed  and  is  documented  in  Reference  [5].  The  reader  is  referred  to 
these  references  for  a  thorough  discussion  of  the  theory  and  practice  of  PSE  methodology. 


4  Test  Problem 


As  implied  in  the  introduction,  the  computational  experiment  is  an  approximate  analog  to 
the  quiet  wind-tunnel  experiment  of  Cattafesta  et  al.  [3].  In  this  section,  we  compare  and 
contrast  the  physical  and  numerical  experiments. 


4.1  The  physical  experiment 

The  reader  is  referred  to  Cattafesta  et  al.  [3],  which  is  summarized  briefly  here.  A  15- 
inch-long  model  of  a  wing  section  with  a  leading-edge  sweep  angle  =  77.1^  is  being 
tested  in  NASA  Langley’s  Mach  3.5  quiet  wind-tunnel  to  investigate  crossflow  instability 
and  transition.  A  three-dimensional  view  of  the  model  is  shown  in  Fig.  2  of  Cattafesta  et  al. 
[3].  The  model  was  designed  originally  to  experience  flow  similar  to  that  on  the  70-degree 
swept  wing  of  an  F-16XL,  which  is  undergoing  flight  experiments  by  N.A.SA.  The  model  cross 
section  is  geometrically  similar  to  that  of  the  first  6.25  percent  chord  of  the  wing  glove  on 
the  modified  F-16XL.  However,  the  sweep  angle  on  the  wind-tunnel  model  was  increased 
to  77  degrees  to  match  the  leading-edge-normal  Mach  number  of  the  Mach  3.5  wind-tunnel 
experiment  to  that  of  the  Mach  2.4  flight  experiment,  for  which  the  normal  Mach  number  is 
0.78. 

Because  crossflow-dominated  transition  is  quite  sensitive  to  controlled  and  random  influ¬ 
ences,  it  was  deemed  necessary  to  conduct  the  experiment  in  a  quiet  facility.  In  the  physical 
experiment,  the  freestream  unit  Reynolds  numbers  were  varied  from  1.5  to  8.0  million  per 
foot  by  varying  the  tunnel  stagnation  pressure,  and  the  angle  of  attack  was  incremented 
from  -2  to  5  degrees.  In  support  of  the  experiment,  Euler  and  Navier-Stokes  calculations 
were  conducted,  as  were  A^-factor  studies  using  the  stability  code  COS.A.L  (Iyer  et  al.  [15]). 
As  in  most  theoretical  studies  of  crossflow  instability,  traveling  waves  were  predicted  to  have 
the  highest  growth  rates;  peak  A'-factors  were  observed  for  frequencies  in  the  range  of  40 


8 


to  60kHz.  The  transition  front  on  the  model  was  estimated  on  the  basis  of  recovery  fac¬ 
tors  obtained  from  surface  thermocouples.  More  recently,  temperature-sensitive  paint  has 
been  used  to  more  finely  resolve  the  transition  front.  The  measured  and  predicted  transition 
fronts  correlated  approximately  for  N  =  14.  However,  regions  on  the  model  with  microscopic 
surface  scratches  showed  the  telltale  streaks  of  stationary  crossflow  vortices. 


4.2  The  numerical  experiment 

Data  for  the  numerical  experiment  were  inferred  from  an  Euler  calculation  of  the  flow  on 
the  wing  model  by  Iyer  [14],  for  which  the  freestream  conditions  were 
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Fig.  6  shows  isobars  on  the  upper  and  lower  surfaces  of  the  model  for  the  freestream 
conditions  given  above.  Near  the  aft  stations  on  the  wing,  the  wing  sections  are  similar 
(Cattafesta  et  al.  [3]),  and  the  isobars  are  nearly  parallel,  which  suggests  that  the  flow  can 
be  approximated  as  quasi-three-dimensional  (although  Cattafesta  et  al.  caution  against  this). 
Accordingly,  we  consider  a  wing  section  perpendicular  to  the  leading  edge  at  a  station  1.091 
ft.  along  the  leading  edge,  as  shown  in  Fig.  7.  The  chord  length  c*  normal  to  the  leading 
edge  at  the  section  of  interest  is  0.2496  ft.  The  surface  pressure  coefficients  interpolated 
from  the  Euler  grid  to  the  wing  section  are  shown  in  Fig.  8.  The  surface  pressure  coefficients 
were  then  used  as  input  data  for  the  boundary-layer  code  WINGBL2  (described  earlier) 
to  derive  the  base  state.  Some  additional  interpretation  of  the  Euler  data,  however,  was 
necessary  to  provide  post-shock  reference  conditions  to  the  boundary- layer  code.  From  the 
Euler  data,  the  post-shock  Mach  number  at  the  point  of  maximum  wingspan  was  taken  to 
be  3.29.  Because  of  the  outward  turning  of  the  flow  through  the  shock,  the  effective  wing 
sweep  angle  increased  by  a  few  degrees  over  the  physical  value.  Finally,  the  post-shock 
reference  temperature  was  adjusted  to  force  agreement  between  the  Euler  and  boundary-layer 
stagnation  temperatures.  These  self-consistent  post-shock  reference  values  are  summarized 
as  follows: 
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=  80.81° 


From  the  geometry,  reference  values,  and  pressure  coefficients,  the  boundary- layer  solu¬ 
tion  was  computed  based  on  the  approximation  of  isentropic  post-shock  flow.  The  boundary- 
layer  and  Euler  solutions  for  the  edge  Mach  number  and  edge  temperature  are  compared 
in  Figs.  9  and  10,  respectively.  The  small  disagreement  between  the  two  solutions  can  be 
attributed  to  two  sources.  First,  although  the  shock  is  weak,  shock  strength  varies  in  az¬ 
imuth  due  to  the  asymmetry  of  the  body.  Consequently,  the  post-shock  flow  is  not  strictly 
isentropic.  Second,  of  course,  the  base  flow  is  not  perfectly  two-dimensional.  Nevertheless, 
the  close  agreement  between  the  two  solutions  suggests  that  the  simplifying  assumptions  are 
reasonable. 

Growth  of  the  boundary-layer  displacement  thickness,  obtained  from  WINGBL2,  is  pre¬ 
sented  in  Fig.  11.  The  favorable  pressure  gradient  induces  a  boundary  layer  that  grows 
linearly,  except  during  the  most  rapid  acceleration  near  the  leading  edge.  The  boundary- 
layer  solution  is  compared  at  two  stations,  .r^  =  l(^*n  —  9-0  =  0.173.  in  Figs.  12 

and  13,  which  show  velocity  and  temperature  profiles,  respectively.  By  definition,  of  course, 
the  streamwise  velocity  vanishes  along  the  attachment  line,  which  corresponds  (nominally) 
to  Xc  =  0.0.  Figure  14  shows  the  solution  also  at  Xc  =  0.173,  but  in  a  coordinate  system 
oriented  along  the  inviscid  streamline.  In  the  rotated  coordinated  system,  u*  and  u*  de¬ 
note  the  tangential  and  crossflow  velocities,  respectively.  Note  the  inflectional  nature  of  the 
crossflow  velocity  (the  component  perpendicular  to  the  inviscid  streamline).  The  crossflow 
Reynolds  number,  defined  in  Reed  and  Saric  [28],  is  shown  as  a  function  of  Xc  in  Fig.  1-5. 
The  maximum  crossflow  Reynolds  number  is  attained  at  approximately  20  percent  chord 
and  decreases  gradually  thereafter. 

Finally,  the  boundary-layer  solution  is  interpolated  onto  the  DNS  grid  by  the  interpolation 
code  INTBL  described  previously.  The  computational  domain  for  the  DNS  spans  approxi¬ 
mately  from  1  to  70  percent  chord  in  streamwise  extent  and  from  the  wall  to  s  =  z*  jS*  —  25 
in  wall-normal  extent.  At  the  inflow  boundary,  6*  =  0.00053  ft.  .A.n  algebraic  mapping 
is  used  in  the  wall-normal  direction  to  cluster  points  close  to  the  wall.  A  second  (lin¬ 
ear)  mapping  is  used  to  remove  the  major  effects  of  the  boundary-layer  growth  shown  in 
Fig.  11,  so  that  the  boundary  layer  remains  of  nearly  constant  thickness  in  the  scaled  coor¬ 
dinate  z*/^*(x*).  iV-factor  studies  by  Chang  [8].  who  used  PSE  methodology,  determined 
that  a  stationary  crossflow  mode  of  spanwise  wavelength  A“  =  10  mm  (0.0328  ft.)  was 
strongly  amplified.  Accordingly,  the  fundamental  spanwise  wavenumber  was  chosen  to  be 
=  0*8*  =  2'k8*  jXy  =  0.1018.  The  width  of  the  computational  domain  was  one  spanwise 
wavelength;  that  is,  1  x  A*  =  0.0328  ft.  The  suction  strip  spanned  from  6.7  to  9.6  percent 
chord,  and  rc^vall  ~  10““*. 


10 


5  Results 


The  DNS  calculation  was  made  with  a  resolution  of  577  x  4  x  97  in  the  streamwise,  spanwise, 
and  wall-normal  directions,  respectively.  The  streamwise  and  wall-normal  resolutions  were 
based  on  the  DNS  experience  of  Joslin  [16]  and  of  Pruett  and  Chang  [27].  The  spanwise 
resolution  was  sufficient  to  resolve  only  the  mean  and  a  single  crossflow  mode;  hence,  dealias¬ 
ing  was  exploited  (in  Fourier  space)  in  the  spanwise  direction  to  remove  energy  in  the  first 
harmonic  of  the  spanwise  fundamental. 

To  establish  the  crossflow  disturbances  throughout  the  domain  required  an  integration  in 
time  to  approximately  t  =  1800,  at  which  time  the  leading  wavefront  had  encountered  the 
outflow  boundary.  Because  the  present  outflow  conditions  are  not  completely  satisfactory,  a 
fraction  of  the  incident  energy  was  reflected  by  the  boundary  and  contaminated  the  upstream 
solution.  Hence,  the  computation  was  halted  prior  to  the  exit  of  the  leading  wavefront  from 
the  domain.  Fig.  16  shows  the  evolution  of  the  maximum  of  the  (dealiased)  perturbation 
spanwise  velocity  as  a  function  of  the  streamwise  coordinate.  Superimposed  on  the  plot  is 
the  amplitude  of  the  r  component  of  the  crossflow  instability  as  predicted  by  linear  PSE 
methodology.  The  initial  amplitude  of  the  PSE  result  is  arbitrarily  scaled  for  the  purpose 
of  comparison.  In  the  vicinity  of  the  suction  strip,  where  receptivity  is  enhanced,  spatial 
transients  are  observed.  Downstream  of  the  suction  strip,  apparently  a  stationary  crossflow 
mode  is  established,  which  grows  in  close  agreement  with  the  prediction  of  linear  PSE.  The 
leading  wave  packet  is  broad  and  spans  approximately  0.35  <  Xc  <  0.7  at  termination  of  the 
calculation.  Moreover,  by  comparison  with  the  linear  PSE  result,  it  appears  that  the  leading 
wave  packet  experiences  amplitudes  one  to  two  orders  of  magnitude  larger  than  the  crossflow 
instability  that  is  deposited  by  its  passing.  Thus,  despite  the  initially  low  amplitude  of  the 
disturbance,  strong  nonlinearities  are  encountered  as  the  leading  w'avefront  passes.  Indeed, 
in  a  previous  simulation  with  a  different  buffer-domain  treatment,  in  which  the  wavefront 
completely  exited  the  computational  domain,  nonlinearities  and  their  associated  Lighthill 
stresses  were  strong  enough  in  the  wave-packet  region  to  generate  significant  acoustic  energy, 
as  shown  in  Fig.  17.  Following  passage  of  the  wavefront,  the  acoustic  radiation  vanished. 

Figure  16  suggests  that  a  linear  crossflow  mode  has  been  established  in  the  DNS  calcu¬ 
lation  in  the  region  0.10  <  x^  <  0.35.  Further  confirmation  of  this  is  evidenced  in  Fig.  18. 
which  compares  the  amplitude  of  the  fundamental  spanwise  Fourier  harmonic  of  the  DNS 
calculation  with  the  amplitude  of  the  crossflow  mode  predicted  by  the  PSE  method  at  the 
station  x^  =  0.22.  For  the  comparison,  the  PSE  and  DNS  results  are  each  normalized  so  that 
the  maximum  spanwise  perturbation  velocity  is  unity.  The  predictions  of  the  two  methods 
agree  very  well  for  all  perturbation  components.  Consequently  for  the  DNS  calculation,  the 
principal  effect  of  the  receptivity  of  the  flow  to  the  simulated  roughness  element  is  the  gen¬ 
eration  of  a  linear  crossflow  mode  immediately  downstream  of  the  suction  strip,  an  expected 
result. 

Fig.  19  depicts  the  alignment  of  the  crossflow  vortices  immediately  downstream  of  the 
suction  strip  in  the  DNS  calculation.  The  disturbances  are  made  visible  as  contours  of  con- 
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stant  disturbance  density  in  a  plane  approximately  1.5  displacement  thicknesses  from  the 
wall.  By  slicing  Fig.  19  in  a  specific  spanwise  plane,  the  vortex  alignment  angle  can  be 
estimated  on  the  basis  of  the  ratio  of  the  wavelengths  in  the  spanwise  and  the  streamwise 
directions.  Whereas  the  dimensional  spanwise  wavelength  is  fixed,  the  streamwise  wave¬ 
length  A*  varies,  as  shown  in  Fig.  20.  Relative  to  the  x  axis,  the  vortex  alignment  angle 
is  tan“^(a*//?*),  where  the  streamwise  wavenumber  a*  =  2t:IXI.  Shown  in  Fig.  21  are  the 
vortex  alignment  angles  as  computed  by  linear  stability  theory  (which  assumes  the  flow  to 
be  locally  parallel)  and  by  the  (nonparallel)  PSE  method.  For  comparison,  the  angle  of  the 
local  inviscid  streamline  relative  to  the  ::  axis  is  also  shown.  The  inviscid  streamline  angle 
is  derived  from  the  boundary-layer  solution  as  tan~^(ug/i<g).  The  vortex  angle  at  Xc  =  0.2 
as  derived  from  the  DNS  results  is  also  shown  and  agrees  closely  with  the  PSE  result.  Both 
the  DNS  and  PSE  approaches  clearly  show  the  tendency  of  crossflow  vortices  to  align  nearly 
along  inviscid  streamlines.  The  maximum  deviation  in  the  streamline  and  vortex  alignment 
angles  is  approximately  three  degrees  at  Xc  =  0.13. 


6  Conclusions 

•  Direct  numerical  simulation  of  crossflow  instability  in  compressible  flows  on  swept 
wings  is  computationally  intensive.  The  simulation  discussed  in  the  previous  section 
required  100  hours  of  CPU  time  on  a  Cray  Y-MP,  despite  relatively  coarse  streamwise 
and  spanwise  resolution.  Several  factors  contribute  to  this  expense.  First,  long-time 
integration  is  required  to  established  the  instability.  Whereas  the  use  of  implicit  or 
semi-implicit  methods  of  time  advancement  could  reduce  the  computational  expense  by 
allowing  longer  time  steps,  such  savings  are  not  guaranteed.  Joslin  [16]  implies  that  the 
time  step  needed  to  suppress  temporal  transients  and  numerical  error  is  considerably 
less  than  that  permitted  by  the  stability  of  his  semi-implicit  approach.  Second,  because 
even  small  disturbances  quickly  generate  a  large-amplitude  leading  wave  packet,  fine 
grid  resolution  is  needed  to  resolve  nonlinearities,  even  if  the  goal  is  to  investigate 
linear  disturbances. 

•  At  the  present  stage  of  development,  additional  effort  is  needed  to  refine  both  the  inflow 
and  outflow  boundary  conditions.  It  is  not  presently  known  whether  the  reflections 
experienced  at  the  outflow  boundary  are  physical  or  numerical  in  origin;  however,  given 
our  previous  experience  with  boundary  difficulties  (Pruett  et  al.  [26]),  we  believe  the 
problem  can  be  rectified  given  sufficient  attention. 

•  Despite  the  limitations  of  the  numerical  experiment,  the  qualitative  and  quantitative 
agreement  between  the  DNS  and  PSE  results  is  good.  The  DNS  confirms  that  the 
crossflow  modes  predicted  to  be  unstable  by  PSE  are  indeed  unstable.  Moreover,  the 
rate  of  growth,  the  the  vortex  orientation  angles,  and  the  modal  structures  predicted  by 
the  two  methods  are  shown  to  agree  closely.  Although  more  validation  is  warranted, 
particularly  of  strongly  nonlinear  interactions,  we  believe  that,  given  the  agreement 
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between  the  PSE  and  DNS  methods,  PSE  methodology  can  be  used  reliably  and  inex¬ 
pensively  to  investigate  stationary  crossflow  instabilities  on  swept-wing  configurations 
relevant  to  the  design  of  the  HSCT. 

•  An  unexpected  observation  was  the  radiation  of  acoustic  energy  from  the  (highly  non¬ 
linear)  developing  wavefront  as  the  crossflow  instability  was  being  established.  Al¬ 
though  this  phenomenon  is  likely  an  artifact  of  the  method  used  to  generate  distur¬ 
bances,  it  underscores  the  capability  of  the  DNS  algorithm  and  suggests  that  the 
method  could  readily  be  adapted  to  applications  in  the  area  of  aeroacoustics. 
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Figures 


.  Individual  terms  of  continuity  ecpation  and  total  residual  at  attach¬ 
ment  line. 
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2.  Individual  terms  of  continuity  equation  and  total  residual  at  x*  =  0.04 
ft. 
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wall-normal  index, 
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5.  Wall-normal  index  of  global  maximum  residual  versus  time. 
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77-Degree  Swept  Wing  in  Mach  3.5  Flow 
Surface  Pressure  Coefficient 


Lower  Surface 


Upper  Surface 


6.  Euler  solution  displaying  isobars  of  constant  pressure  on  upper  and 
lower  surfaces  of  wing  model  of  Cattafesta  et  al,  [3] 
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Wall-normal  distributions  of  streamwise  and  spanwise  velocities  at  at¬ 
tachment  line  and  at  Xc  =  0.173  on  upper  wing  surface. 


temperature  (Rankine) 
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13.  Wall-normal  distribution  of  temperature  at  attachment  line  and  at  Xr  = 

0.173. 
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77-Degree  Swept  Wing 
Mach  3.29 


Disurbance  Density 


Acoustic  Radiation 
from  Developing  Wavefront 


Contours  of  constant  disturbance  density  in  computational  plane  nor¬ 
mal  to  both  leading  edge  and  wing  surface  showing  acoustic  energy 
generated  by  strong  nonlinearity  in  wave-packet  region. 


amplitude 
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18.  Comparison  of  amplitude  distributions  of  components  of  fundamen¬ 
tal  spanwise  Fourier  component  of  DNS  calculation  with  PSE-derived 
crossflow  disturbance  at  Xc  =  0.22.  In  each  case,  results  are  normal¬ 
ized  so  that  maximum  spanwise  perturbation  velocity  is  unity.  (Even 
indexed  DNS  values  omitted.) 
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19.  Contours  of  constant  perturbation  density  in  surface  approximately  1.5 
displacement  thicknesses  from  wall  showing  alignment  of  crossflow  vor¬ 
tices.  Computational  domain  is  replicated  once  in  span  for  aesthetics. 
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21.  Comparison  of  vortex-alignment  angle  and  inviscid-streamline  angle 
relative  to  rc-axis  showing  orientation  of  crossflow  mode  in  direction  of 
inviscid  streamlines. 
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